{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Q4qew2DscIPu"
   },
   "source": [
    "# 第19章 马尔可夫链蒙特卡罗法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 蒙特卡罗法是通过基于概率模型的抽样进行数值近似计算的方法，蒙特卡罗法可以用于概率分布的抽样、概率分布数学期望的估计、定积分的近似计算。\n",
    "\n",
    "随机抽样是蒙特卡罗法的一种应用，有直接抽样法、接受拒绝抽样法等。接受拒绝法的基本想法是，找一个容易抽样的建议分布，其密度函数的数倍大于等于想要抽样的概率分布的密度函数。按照建议分布随机抽样得到样本，再按要抽样的概率分布与建议分布的倍数的比例随机决定接受或拒绝该样本，循环执行以上过程。\n",
    "\n",
    "马尔可夫链蒙特卡罗法数学期望估计是蒙特卡罗法的另一种应用，按照概率分布$p(x)$抽取随机变量$x$的$n$个独立样本，根据大数定律可知，当样本容量增大时，函数的样本均值以概率1收敛于函数的数学期望\n",
    "\n",
    "$$\\hat { f } _ { n } \\rightarrow E _ { p ( x ) } [ {f ( x )} ] , \\quad n \\rightarrow \\infty$$\n",
    "\n",
    "\n",
    "计算样本均值$\\hat { f } _ { n } $，作为数学期望$E _ { p ( x ) } [ {f ( x )} ] $的估计值。\n",
    "\n",
    "2. 马尔可夫链是具有马尔可夫性的随机过程$$P ( X _ { t } | X _ { 0 } X _ { 1 } \\cdots X _ { t - 1 } ) = P ( X _ { t } | X _ { t - 1 } ) , \\quad t = 1,2 , \\cdots$$\n",
    "通常考虑时间齐次马尔可夫链。有离散状态马尔可夫链和连续状态马尔可夫链，分别由概率转移矩阵$P$和概率转移核$p(x,y)$定义。\n",
    "\n",
    "满足$\\pi = P \\pi$或$\\pi ( y ) = \\int p ( x , y ) \\pi ( x ) d x$的状态分布称为马尔可夫链的平稳分布。\n",
    "\n",
    "马尔可夫链有不可约性、非周期性、正常返等性质。一个马尔可夫链若是不可约、非周期、正常返的，则该马尔可夫链满足遍历定理。当时间趋于无穷时，马尔可夫链的状态分布趋近于平稳分布，函数的样本平均依概率收敛于该函数的数学期望。\n",
    "\n",
    "$$\\operatorname { lim } _ { t \\rightarrow \\infty } P ( X _ { t } = i | X _ { 0 } = j ) = \\pi _ { i } , \\quad i = 1,2 , \\cdots ; \\quad j = 1,2$$\n",
    "$$\\hat { f } _ { t } \\rightarrow E _ { \\pi } [ {f ( X )} ] , \\quad t \\rightarrow \\infty$$\n",
    "可逆马尔可夫链是满足遍历定理的充分条件。\n",
    "\n",
    "3. 马尔可夫链蒙特卡罗法是以马尔可夫链为概率模型的蒙特卡罗积分方法，其基本想法如下：\n",
    "\n",
    "（1）在随机变量$x$的状态空间$\\chi$上构造一个满足遍历定理条件的马尔可夫链，其平稳分布为目标分布$p(x)$；\n",
    "\n",
    "（2）由状态空间的某一点$X_0$出发，用所构造的马尔可夫链进行随机游走，产生样本序列$X _ { 1 } , X _ { 2 } , \\cdots , X _ { t } , \\cdots$；\n",
    "\n",
    "（3）应用马尔可夫链遍历定理，确定正整数$m$和$n(m<n)$，得到样本集合$\\{ x _ { m + 1 } , x _ { m + 2 } , \\cdots , x _ { n } \\}$，进行函数$f(x)$的均值（遍历均值）估计：\n",
    "$$\\hat { E } f = \\frac { 1 } { n - m } \\sum _ { i = m + 1 } ^ { n } {f ( x _ { i } )}$$\n",
    "\n",
    "4. Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法。假设目标是对概率分布$p ( x )$进行抽样，构造建议分布$q ( x , x ^ { \\prime } )$，定义接受分布$\\alpha ( x , x ^ { \\prime } )$进行随机游走，假设当前处于状态$x$，按照建议分布$q ( x , x ^ { \\prime } )$机抽样，按照概率$\\alpha ( x , x ^ { \\prime } )$接受抽样，转移到状态 $x ^ { \\prime }$，按照概率$1- \\alpha ( x , x ^ { \\prime } )$拒绝抽样，停留在状态$x$，持续以上操作，得到一系列样本。这样的随机游走是根据转移核为$p ( x , x ^ { \\prime } ) = q ( x , x ^ { \\prime } ) \\alpha ( x , x ^ { \\prime } )$的可逆马尔可夫链（满足遍历定理条件）进行的，其平稳分布就是要抽样的目标分布$p ( x )$。\n",
    "\n",
    "5. 吉布斯抽样（Gibbs sampling）用于多元联合分布的抽样和估计吉布斯抽样是单分量 Metropolis-Hastings-算法的特殊情况。这时建议分布为满条件概率分布$$q ( x , x ^ { \\prime } ) = p ( x _ { j } ^ { \\prime } | x _ { - j } )$$\n",
    "吉布斯抽样的基本做法是，从联合分布定义满条件概率分布，依次从满条件概率分布进行抽样，得到联合分布的随机样本。假设多元联合概率分布为$p ( x ) = p ( x _ { 1 } , x _ { 2 } , \\cdots , x _ { k } )$，吉布斯抽样从一个初始样本$x ^ { ( 0 ) } = ( x _ { 1 } ^ { ( 0 ) } , x _ { 2 } ^ { ( 0 ) } , \\cdots , x _ { k } ^ { ( 0 ) } ) ^ { T }$出发，不断进行迭代，每一次迭代得到联合分布的一个样本$x ^ { ( i ) } = ( x _ { 1 } ^ { ( i ) } , x _ { 2 } ^ { ( i ) } , \\cdots , x _ { k } ^ { ( i ) } ) ^ { T }$，在第$i$次迭代中，依次对第$j$个变量按照满条件概率分布随机抽样$p ( x _ { j } | x _ { 1 } ^ { ( i ) } , \\cdots , x _ { j - 1 } ^ { ( i ) },x _ { j + 1 } ^ { ( i - 1 ) } , \\cdots , x _ { k } ^ { ( i - 1 ) } ) , j = 1,2 , \\cdots , k$，得到$x _ { j } ^ { ( i ) }$最终得到样本序列$\\{ x ^ { ( 0 ) } , x ^ { ( 1 ) } , \\cdots , x ^ { ( n ) } \\}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qKLGsUc-cNIM"
   },
   "source": [
    "---\n",
    "蒙特卡洛法（Monte Carlo method) ， 也称为统计模拟方法 (statistical simulation method) ， 是通过从概率模型的随机抽样进行近似数值计\n",
    "\n",
    "算的方法。 马尔可夫链陟特卡罗法 (Markov Chain Monte Carlo, MCMC)， 则是以马尔可夫链 (Markov chain）为概率模型的蒙特卡洛法。\n",
    "\n",
    "马尔可夫链蒙特卡罗法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布， 首先基于该马尔可夫链进行随机游走， 产生样本的序列，\n",
    "\n",
    "之后使用该平稳分布的样本进行近似数值计算。\n",
    "\n",
    "Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法，Metropolis等人在 1953年提出原始的算法，Hastings在1970年对之加以推广，\n",
    "\n",
    "形成了现在的形式。吉布斯抽样(Gibbs sampling)是更简单、使用更广泛的马尔可夫链蒙特卡罗法，1984 年由S. Geman和D. Geman提出。\n",
    "\n",
    "马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习\n",
    "\n",
    "与推理，是重要的统计学习计算方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jE5JCrmed0C4"
   },
   "source": [
    "一般的蒙特卡罗法有**直接抽样法**、**接受-拒绝抽样法**、 **重要性抽样法**等。\n",
    "\n",
    "接受-拒绝抽样法、重要性抽样法适合于概率密度函数复杂 （如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂），不能直接抽样的情况。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2ZBr1H0Md0gv"
   },
   "source": [
    "### 19.1.2 数学期望估计\n",
    "\n",
    "一舣的蒙特卡罗法， 如直接抽样法、接受·拒绝抽样法、重要性抽样法， 也可以用于数学期望估计 （estimation Of mathematical expectation)。\n",
    "\n",
    "假设有随机变量$x$， 取值 $x\\in X$, 其概率密度函数为 $p(x)$, $f(x)$ 为定义在 $X$ 上的函数， 目标是求函数 $f(x)$ 关于密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)]$。\n",
    "\n",
    "\n",
    "针对这个问题，蒙特卡罗法按照概率分布 $p(x)$ 独立地抽取 $n$ 个样本$x_{1}, x_{2},...,x_{n}$，比如用以上的抽样方法，之后计算函\n",
    "\n",
    "数$f(x)$的样本均值$\\hat f_{n}$  \n",
    "\n",
    "$\\hat f_{n} = \\frac{1} {n}\\sum_{i=1}^{n}f(x_{i})$\n",
    "\n",
    "\n",
    "作为数学期望$E_{p(x)}[f(x)]$近似值。\n",
    "\n",
    "根据大数定律可知， 当样本容量增大时， 样本均值以概率1收敛于数学期望：  \n",
    "\n",
    "$\\hat f_{n} \\rightarrow E_{p(x)}[f(x)], n \\rightarrow \\infty $\n",
    "\n",
    "这样就得到了数学期望的近似计算方法：  \n",
    "\n",
    "$E_{p(x)}[f(x)] \\approx \\frac{1} {n}\\sum_{i=1}^{n}f(x_{i})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZtV3LYrUh9jG"
   },
   "source": [
    "### 马尔可夫链 \n",
    "\n",
    "考虑一个随机变量的序列 $X = {X_{0}, X_{1},..., X(t),...}$ 这里 $X_{t}$，表示时刻 $t$ 的随机变量， $t = 0, 1, 2...$. \n",
    "\n",
    "每个随机变量 $X_{t}(t=0,1,2,...)$ 的取值集合相同， 称为状态空间， 表示为$S$.  随机变量可以是离散的， 也可以是连续的。\n",
    "\n",
    "以上随机变量的序列构成随机过程（stochastic process)。\n",
    "\n",
    "假设在时刻 $0$ 的随机变量 $X_{0}$ 遵循概率分布 $P(X_{0}) = \\pi$，称为初始状态分布。在某个时刻 $t>=1$ 的随机变量 $X_{t}$与前\n",
    "\n",
    "一个时刻的随机变量 $X_{t-1}$ 之间有条件分布 $P(X_{t}|X_{t-1})$ 如果 $X_{t}$ 只依赖于 $X_{t-1}$, 而不依赖于过去的随机变量 \n",
    "\n",
    "${X_{0}，X_{1},...，X_{t-2}}$ 这一性质称为马尔可夫性，即  \n",
    "\n",
    "$P(X_{t}|X_{0},X_{1},...,X_{t-1}) = P(X_{t}|X_{t-1}), t=1,2,...$\n",
    "\n",
    "具有马尔可夫性的随机序列$X = {X_{0}, X_{1},..., X(t),...}$称为马尔可夫链， 或马尔可夫过程（Markov process)。 条件概率分布 \n",
    "\n",
    "$P(X_{t}|X_{t-1})$ 称为**马尔可夫链的转移概率分布**。 **转移概率分布决定了马尔可夫裢的特性**。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BL8oLbgBttgs"
   },
   "source": [
    "#### 平稳分布  \n",
    "\n",
    "设有马尔可夫链$X = {X_{0}, X_{1},..., X(t),...}$，其状态空间为 $S$,转移概率矩阵为 $P=(p_{ij})$， 如果存在状态空间 $S$ 上的一个分布  \n",
    "\n",
    "$\\pi = \\begin{bmatrix}\n",
    "\\pi_{1}\\\\ \n",
    "\\pi_{2}\\\\ \n",
    "\\vdots \\end{bmatrix}$\n",
    "\n",
    "使得  \n",
    "\n",
    "$\\pi = P\\pi$\n",
    "\n",
    "则称丌为马尔可夫裢$X = {X_{0}, X_{1},..., X(t),...}$的平稳分布。\n",
    "\n",
    "\n",
    "直观上，如果马尔可夫链的平稳分布存在，那么以该平稳分布作为初始分布，面向未来进行随机状态转移，之后任何一个时刻的状态分布都是该平稳分布。\n",
    "\n",
    "**引理19．1**\n",
    "\n",
    "给定一个马尔可夫链$X = {X_{0}, X_{1},..., X(t),...}$, 状态空间为$S$, 移概率矩阵为$P=(p_{ij})$， 则分布 $\\pi=(\\pi_{1}, \\pi_{2},...)^{T}$ 为 $X$ 的平稳分布的充要条件是$\\pi=(\\pi_{1}, \\pi_{2},...)^{T}$是下列方程组的解：\n",
    "\n",
    "$x_{i} = \\sum_{j}p_{ij}x_{j}, i=1,2,...$  \n",
    "\n",
    "$x_{i} >= 0, i = 1,2,...$  \n",
    "\n",
    "$\\sum_{i}x_{i} = 1$  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HqSHX7PSwOjP"
   },
   "source": [
    "### 吉布斯采样\n",
    "\n",
    "输入： 目标概率分布的密度函数$p(x)$, 函数$f(x)$;\n",
    "\n",
    "输出： $p(x)$的随机样本 $x_{m+1}, x_{m+2}, ..., x_{n}$，函数样本均值 $f_{mn}$;\n",
    "\n",
    "参数： 收敛步数$m$, 迭代步数 $n$.\n",
    "\n",
    "\n",
    "1. 初始化。给出初始样本 $x^{0} = $($x^{0}_{1}, x^{0}_{2},..., x^{0}_{k}$)$^{T}$.\n",
    "\n",
    "2. 对$i$循环执行  \n",
    " 设第$i-1$次迭代结束前的样本为$x^{i-1} = $($x^{i-1}_{1}, x^{i-1}_{2},..., x^{i-1}_{k}$)$^{T}$，则第$i$次迭代进行如下几步操作：  \n",
    "\n",
    " + (1)由满条件分布 $p(x_{1}|x^{i-1}_{2},...,x^{i-1}_{k})$ 抽取 $x^{i}_{1}$  \n",
    " \n",
    " + ...\n",
    " \n",
    " + (j)由满条件分布 $p(x_{j}|x^{i}_{1},...,x^{i}_{j-1}, x^{i-1}_{j+1},..., x^{i-1}_{k})$ 抽取 $x^{i}_{j}$   \n",
    " \n",
    " + (k)由满条件分布 $p(x_{k}|x^{i}_{1},...,x^{i}_{k})$ 抽取 $x^{i}_{k}$ \n",
    " \n",
    "得到第 $i$ 次迭代值 $x^{(i)} = (x^{(i)}_{1}, x^{(i)}_{2},..., x^{(i)}_{k})^{T}$.\n",
    "\n",
    "\n",
    " 3. 得到样本集合\n",
    "    \n",
    " {$x^{(m+1)}, x^{(m+2)},..., x^{(n)}$}\n",
    " \n",
    " 4. 计算\n",
    " \n",
    " $f_{mn} = \\frac{1}{n-m}\\sum_{i=m+1}^{n}f(x^{(i)})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5ZeiXcVWBQZb"
   },
   "source": [
    "--------------------------------------------------------------------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PIcMXLgQBSio"
   },
   "source": [
    "#### 网络资源：\n",
    "\n",
    "LDA-math-MCMC 和 Gibbs Sampling： https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling  \n",
    "\n",
    "MCMC蒙特卡罗方法： https://www.cnblogs.com/pinard/p/6625739.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "kIIlKmr0I8d_",
    "outputId": "265030d7-7e28-443c-eed0-e63ee8b21a73"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.23076935 0.30769244 0.46153864]]\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "\n",
    "transfer_matrix = np.array([[0.6, 0.2, 0.2], [0.3, 0.4, 0.3], [0, 0.3, 0.7]],\n",
    "                           dtype='float32')\n",
    "start_matrix = np.array([[0.5, 0.3, 0.2]], dtype='float32')\n",
    "\n",
    "value1 = []\n",
    "value2 = []\n",
    "value3 = []\n",
    "for i in range(30):\n",
    "    start_matrix = np.dot(start_matrix, transfer_matrix)\n",
    "    value1.append(start_matrix[0][0])\n",
    "    value2.append(start_matrix[0][1])\n",
    "    value3.append(start_matrix[0][2])\n",
    "print(start_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 269
    },
    "colab_type": "code",
    "id": "A2oybaKXJGqd",
    "outputId": "3b2abcf8-71f1-4953-8536-667de8a1f36e"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8ddnJpM9ELIiEPYACQHBBlGxIqK4At5W/bmWan9X/VWrrT8XvFq1tP706m3rbS9X6+21thc3ipWiaFUQBGkVEsWFDIEQQcKShIQEsicz398fMwkh6yRMcmb5PB/M4+xnPod58J7D95zzHTHGoJRSKrTZrC5AKaXUwNOwV0qpMKBhr5RSYUDDXimlwoCGvVJKhYEIqwvoKCUlxYwdO9bqMpRSKqjk5+cfMcakdrc84MJ+7Nix5OXlWV2GUkoFFRHZ19NybcZRSqkwoGGvlFJhQMNeKaXCgIa9UkqFAQ17pZQKAxr2SikVBjTslVIqDATcffYq8LmNm2Z3M42uRppcTTS6Gml2NeMyLlrcLbSYFlxu14lpdwsu48LldtFiWnAbN27jxhjjGccz7jKuTvPcxo3BYIzpcQicGLab16q1K+/267Qt67C9d6KTjvvsbX6g0e7MA196XDpXT7p6QPatYR9G6prrKKsro6KhgpqmGo43H6emqYaa5pq24fGm4ydNtw/0tmB3N1t9KKqfBLG6BNWDaanTNOxVz2qbazlUc4jSulLPq9YzPFx3uG38eNPxbrePkAjiI+OJd8STEJlAfGQ8I+JHEGOPIdIeSZQ9ikh75EnjUfYoouxROGwOHHYHEbYIIiQCu9iJsEVgt9k90zY7drHjsDmwiQ2b2LCLHZvYEBHPPDzjdrG3zRMEEWkbtq4DnDTdPsDalnu3aR1vPzwxaLcdctL2HZd33L9SwUbDPgjVNddReLSQr458xVdHvmJHxQ72Hev8pHRydDLpcelkJGSQm55Lelw66bHpJMckk+DwBHpCZALxjnii7FEaZEqFMA37ANfsambX0V3sqNjhCfeKr9hTtQe3cQOQFptGTnIOiyYsIiMhg/TYdNLj0kmLScNhd1hcvVIqUGjYB6BmVzMflnzI67tf55NDn7S1kSdGJTI1ZSoXZFzA1OSp5KTkkBrbbSd3SinVRsM+gBRXF/PG7jdYs2cNlQ2VpMWmce2Ua5meOp2c5BxGxo/UphalVL9o2FusrrmOd/e+y192/4Xt5duJkAjmZszlO5nfYc6IOdhtdqtLVEqFAJ/CXkQuAf4dsAO/N8Y82c16VwF/BmYZY/JEZCzgBAq9q3xsjLn9VIsOdsYYvjryFa/vfp13vn6HupY6xg4Zyz3fuoeFExaSEpNidYlKqRDTa9iLiB1YDlwElADbRGSNMaagw3oJwF3AJx12sccYM8NP9QY9Z4WTR//+KM5KJzERMSwYs4DvZH6HmWkztYlGKTVgfDmzPxMoMsYUA4jIq8BioKDDej8HngLu9WuFIcIYw8s7X+aXeb9kWPQwfnrWT7ls3GXER8ZbXZpSKgz40jfOSGB/u+kS77w2IjITyDDGvNXF9uNE5DMR+VBEvt3VG4jIrSKSJyJ55eXlvtYeNKoaqrhrw108ufVJ5oyYw+sLX+eayddo0CulBo0vZ/ZdtS20dbIhIjbg18D3u1jvEDDaGFMhIt8CVovIVGPMsZN2ZszzwPMAubm5IdWBR97hPJZuXkplQyVLz1zK9VOu1+YapdSg8yXsS4CMdtOjgIPtphOAHGCjN8SGA2tEZJExJg9oBDDG5IvIHmASEPK/KO5yu3j+y+d57vPnyEjIYMVlK8hOzra6LKVUmPIl7LcBmSIyDjgAXAtc37rQGFMNtN0+IiIbgXu9d+OkApXGGJeIjAcygWI/1h+QSmtLWbp5KXmleSwcv5CHznqIOEec1WUppcJYr2FvjGkRkTuBd/HcevmCMWaHiCwD8owxa3rY/DxgmYi0AC7gdmNMpT8KD1Qf7v+Qh7c8TKOrkcfPfZxFExZZXZJSSiGB1sd1bm6uycsLvlaeJlcTv87/NSucK5iSNIWnz3uasUPHWl2WUipMiEi+MSa3u+X6BK0fGGNYunkp7+97n+unXM89ufcQZY+yuiyllGqjYe8Hq3av4v1973PPt+7h5pybrS5HKaU60d+gPUXFVcU8tfUpzhlxDkumLrG6HKWU6pKG/SlodDVy/6b7iXXE8vi5j2MT/etUSgUmbcY5Bc/kP0Ph0UKWz1+unZcppQKanor20+aSzaxwruCGrBs4b9R5VpejlFI90rDvhyP1R3h4y8NMGjaJn3zrJ1aXo5RSvdJmnD5yGzcPf/Qwtc21vHDxC3qLpVIqKOiZfR+tKFjBloNbuH/W/UxInGB1OUop5RMN+z5wVjj59ae/5oKMC7h60tVWl6OUUj7TsPdRXXMd92+6n6SoJH52zs+0m2KlVFDRNnsfPbXtKfYd28d/LfgvEqMTrS5HKaX6RM/sffDe3vd4fffr3JJzC7NPm211OUop1Wca9r04XHuYx/7xGDnJOdwx8w6ry1FKqX7RsO+B27hZunkpLreLp857CofNYXVJSinVLxr2Pdh6eCv5pfncO+teMoZk9L6BUkoFKA37HqwuWk2CI4GF4xdaXYpSSp0SDftuHG86zvp967l03KVER0RbXY5SSp0SDftuvLv3XRpcDVw58UqrS1FKqVOmYd+N1UWrmTB0AjkpOVaXopRSp0zDvgvF1cV8Xv45iycu1idllVIhQcO+C38t+it2sbNwgl6YVUqFBg37DlxuF2/teYtzR56rvz6llAoZGvYd/P3g3ymrL9MLs0qpkKJh38HqotUkRiUyd9Rcq0tRSim/0bBvp7qxmg37N3D5+Mtx2LVrBKVU6NCwb2dt8Vqa3c3ahKOUCjka9u2sLlrNlKQpTEmaYnUpSinlVxr2XoWVhTgrnXpWr5QKSRr2Xn/d81cibBFcNu4yq0tRSim/07AHmt3NrC1ey7yMeQyLHmZ1OUop5Xca9sCmkk1UNlSyeMJiq0tRSqkBoWGP58JsSkwKc0bOsboUpZQaEGEf9kfqj7C5ZDMLxy8kwhZhdTlKKTUgfAp7EblERApFpEhElvaw3lUiYkQkt928B73bFYrIxf4o2p/WFq/FZVx6F45SKqT1eiorInZgOXARUAJsE5E1xpiCDuslAHcBn7Sblw1cC0wFRgDrRGSSMcblv0PoP2MMq4tWMz1lOuMTx1tdjlJKDRhfzuzPBIqMMcXGmCbgVaCrK5k/B54CGtrNWwy8aoxpNMZ8DRR59xcQCioKKKoqYvFEvTCrlAptvoT9SGB/u+kS77w2IjITyDDGvNXXbb3b3yoieSKSV15e7lPh/vBG0RtE2aO4ZNwlg/aeSillBV/CvqufajJtC0VswK+B/9vXbdtmGPO8MSbXGJObmprqQ0mnrtHVyDtfv8MFoy9gSOSQQXlPpZSyii+3n5QAGe2mRwEH200nADnARu9P+A0H1ojIIh+2tcyG/Rs41nRML8wqpcKCL2f224BMERknIpF4LriuaV1ojKk2xqQYY8YaY8YCHwOLjDF53vWuFZEoERkHZAJb/X4U/bC6aDXD44Yze/hsq0tRSqkB12vYG2NagDuBdwEnsNIYs0NElnnP3nvadgewEigA/gbcEQh34pTWlvKPg/9g4fiF2G12q8tRSqkB59NTRMaYt4G3O8x7pJt1z+8w/TjweD/rGxBvFb+F27i1CUcpFTbC8gnavNI8JiZOZPSQ0VaXopRSgyJkwr6pxc32/VWUHmvocT1jDAUVBWQnZw9SZUopZb2QCfuquiauXL6Fv311uMf1yuvLqWyo1LBXSoWVkAn71IQoEmMdFJYe73G9nZU7AchKyhqMspRSKiCETNiLCJPSEth1uOewL6goQBAmJ00epMqUUsp6IRP2AJOGx7Or9DjGdHpIt42zwsmYIWOIc8QNYmVKKWWt0Ar79ASONbRQeqyx23WclU5twlFKhZ2QC3uAXd2021c1VHGo9hBZyRr2SqnwElZh76x0AmjYK6XCTkiFfVJcJCnxUb2HvTbjKKXCTEiFPcDk4fEUltZ0ucxZ4WRE3AiGRg0d5KqUUspaIRf2mWkJ7C49jtvd+Y4cZ6VTm3CUUmEp5MJ+8vAE6ppcHKiqP2l+TVMN+47t0yYcpVRYCrmwn5QeD3S+SFt4tBDQi7NKqfAUcmGf2XZHzsnt9s4KvTirlApfIRf2Q6IdnDY0utOZvbPSSUpMCqmxg/Mbt0opFUh8+vGSYDMpPaHLsNezeqUCT3NzMyUlJTQ09Nw9ufKIjo5m1KhROByOPm0XomEfzz+KK3C5DXab0NDSQHFVMfMy5lldmlKqg5KSEhISEhg7diwiYnU5Ac0YQ0VFBSUlJYwbN65P24ZcMw54zuybWtzsq6gFYPfR3biMS8/slQpADQ0NJCcna9D7QERITk7u1/+CQjLsJw8/+SKtdpOgVGDToPddf/+uQjLsJ6adfPtlQUUBQyKHMCJuhJVlKaWUZUIy7GMjIxidFNsW9jsrd5KVnKVnD0opn33/+99n1apVA7LvzZs3M3XqVGbMmEF9fX23651//vnk5eX55T1DMuzBc5F2V+lxmt3N7Dq6i+wk/c1ZpZT1XC4XL730Evfeey/bt28nJiZmUN43hMM+geLyWgorimh2NzMlaYrVJSmlAtif/vQnpk+fzumnn85NN90EwKZNmzjnnHMYP378SWf5Tz/9NLNmzWL69Ok8+uijbfNXrFjBmWeeyYwZM7jttttwuVwAxMfH88gjjzB79myeeOIJVq5cybJly7jhhhvYuHEjV1xxRds+7rzzTl588UW/H19I3noJnrBvcRs++mY7oBdnlQoGP3tzBwUHj/l1n9kjhvDowqk9rrNjxw4ef/xxtmzZQkpKCpWVldxzzz0cOnSIjz76iJ07d7Jo0SKuuuoq3nvvPXbv3s3WrVsxxrBo0SI2bdpEamoqr732Glu2bMHhcPDDH/6Ql156ie9973vU1taSk5PDsmXLACgqKuKKK67gqquuYuPGjX493u6EdNgDfFq6g9iIWMYMGWNxRUqpQPXBBx9w1VVXkZKSAkBSUhIAV155JTabjezsbEpLSwF47733eO+995g5cyYANTU17N69my+++IL8/HxmzZoFQH19PWlpaQDY7Xa++93vDvZhnSRkw358ahx2m7CnqpApSVOwSci2WCkVMno7Ax8oxpgub+CIioo6aZ3W4YMPPshtt9120rq//e1vWbJkCU888USn/URHR2O327t874iICNxud9v0QD1JHLIJGO2wMzo5miPNX2t7vVKqR/Pnz2flypVUVFQAUFlZ2e26F198MS+88AI1NZ7neA4cOEBZWRnz589n1apVlJWVte1j3759vb73mDFjKCgooLGxkerqatavX++HI+osZM/sAUan1XLE3ajt9UqpHk2dOpWHHnqIuXPnYrfb25pourJgwQKcTidnn3024Ln4umLFCrKzs/nFL37BggULcLvdOBwOli9fzpgxPTchZ2RkcM011zB9+nQyMzN7fO9TIa3/NQkUubm5xl/3ld7x1/9mU9UzvHzpSqalaeArFYicTidZWfrvsy+6+jsTkXxjTG5324RsMw4AkSUYdwTupjSrK1FKKUuFdNgfbfkad+Nwisu7f0JNKaXCQciGvTGGvTW7MI0jOv1qlVJKhZuQDfsDNQc43nScFMd4dh0+3vsGSikVwnwKexG5REQKRaRIRJZ2sfx2EflSRLaLyEciku2dP1ZE6r3zt4vIc/4+gO7srNwJwPihUygs1bBXSoW3XsNeROzAcuBSIBu4rjXM23nZGDPNGDMDeAr4Vbtle4wxM7yv2/1VeG8KKgqwi52Zw6dQcrSe2saWwXprpZQKOL6c2Z8JFBljio0xTcCrwOL2Kxhj2ndmEQdYfj+ns9LJhMQJZA1PBmB3mbbbK6XCly9hPxLY3266xDvvJCJyh4jswXNmf1e7ReNE5DMR+VBEvt3VG4jIrSKSJyJ55eXlfSi/a8YYCioKmJI0hcnprb9apU05Sqnw5UvYd/WLH53O3I0xy40xE4AHgIe9sw8Bo40xM4F7gJdFZEgX2z5vjMk1xuSmpqb6Xn03yuvLqWyoJDs5m4ykWKIibHqRVinVrdraWi6//HJOP/10cnJyeO2111i/fj0zZ85k2rRp3HLLLTQ2Nnba7tChQ5x33nnMmDGDnJwcNm/eDMArr7zCtGnTyMnJ4YEHHhjsw+mSL90llAAZ7aZHAQd7WP9V4FkAY0wj0Ogdz/ee+U8C/POIbDdaL85mJWVhtwmZ6fHs0mYcpQLfO0vh8Jf+3efwaXDpkz2u8re//Y0RI0awdu1aAKqrq8nJyWH9+vVMmjSJ733vezz77LP8+Mc/Pmm7l19+mYsvvpiHHnoIl8tFXV0dBw8e5IEHHiA/P59hw4axYMECVq9ezZVXXunf4+ojX87stwGZIjJORCKBa4E17VcQkcx2k5cDu73zU70XeBGR8UAmUOyPwntSUFGAIExOmgx4ujvWM3ulVHemTZvGunXreOCBB9i8eTN79+5l3LhxTJo0CYAlS5awadOmTtvNmjWLP/zhDzz22GN8+eWXJCQksG3bNs4//3xSU1OJiIjghhtu6HLbwdbrmb0xpkVE7gTeBezAC8aYHSKyDMgzxqwB7hSRC4Fm4CiwxLv5ecAyEWkBXMDtxpjuu5PzE2eFkzFDxhDniAM8Yf+XTw9QXd/M0BjHQL+9Uqq/ejkDHyiTJk0iPz+ft99+mwcffJAFCxZ0ud4nn3zS1rXxsmXL2n64ZO3atdx0003cd999DBnSqaU6IPjU66Ux5m3g7Q7zHmk3fnc3270OvH4qBfaHs9LJjNQZbdOtF2l3lx4nd2zSYJejlApwBw8eJCkpiRtvvJH4+Hiee+459u7dS1FRERMnTuR//ud/mDt3LrNnz2b79u1t2+3bt4+RI0fyz//8z9TW1vLpp5/ywAMPcPfdd3PkyBGGDRvGK6+8wo9+9CMLj84j5Lo4rmqo4lDtIa6bcl3bvMz0eAAKNeyVUl348ssvue+++7DZbDgcDp599lmqq6u5+uqraWlpYdasWdx+e+fHhDZu3MjTTz+Nw+EgPj6eP/3pT5x22mk88cQTzJs3D2MMl112GYsXL+7iXQdXyIW9s9IJnPybsyMTY4iLtLNb+8hRSnXh4osv5uKLL+40/7PPPutxuyVLlrBkyZJO86+//nquv/56v9XnDyHXN05b2CedCHsRITM9gUK9SKuUClOhF/YVTkbEjWBo1NCT5k9OT2B3mYa9Uio8hV7YVzq7/BnCzPR4jtQ0UVHT+cEIpZQKdSEV9jVNNew7tu+kJpxWk4e3dpug7fZKqfATUmFfeLQQoMsze+0jRykVzkIq7J0VnS/OtkpNiGJojEP7tldKhaXQCvtKJykxKaTGdu5MTUQ8F2k17JVSfrJ3715ycnKsLsMnIRf2XZ3Vt8pMj6fw8HGMsby7faWUGlQh81BVQ0sDxVXFzMuY1+06k4cn8NInLZQdbyR9SPQgVqeUCmS1tbVcc801lJSU4HK5+OlPf0phYSFvvvkm9fX1nHPOOfzud79DRMjPz+eWW24hNjaWc8891+rSfRYyYV/TXMOCsQvITc/tdp3MNM9F2sLDxzXslQpA/7r1X9u6KPeXKUlTeODMnvuU76qL44suuohHHvF0AXbTTTfx1ltvsXDhQm6++WZ++9vfMnfuXO677z6/1jqQQqYZJyUmhafOe4qzR5zd7TqTvH3k6B05Sqn2OnZxPHToUDZs2MDs2bOZNm0aH3zwATt27KC6upqqqirmzp0LeL4EgkXInNn7Ijk+ipT4KA17pQJUb2fgA6WrLo6XL19OXl4eGRkZPPbYYzQ0NGCMQaSrH+8LfCFzZu+rSenxFOqDVUqpdg4ePEhsbCw33ngj9957L59++ikAKSkp1NTUsGrVKgASExMZOnQoH330EQAvvfSSZTX3VVid2YPnh0z+nLcft9tgswXnN7RSyr+66uJ49erVTJs2jbFjxzJr1qy2df/whz+0XaDtqqfMQCWBdhtibm6uycsbuJ+offmTb/iXN75k8/3zyEiKHbD3UUr5xul0kpXV/S3TqrOu/s5EJN8Y0+0dKmHXjDN5uOcirfaAqZQKJ2EX9hPbbr/UdnulVPgIu7AfGuPgtKHR2m2CUgEk0JqTA1l//67CLuzBc5FWO0RTKjBER0dTUVGhge8DYwwVFRVER/f9odCwuxsHIGfkEH73YTHHGpoZEu2wuhylwtqoUaMoKSmhvLzc6lKCQnR0NKNGjerzdmEZ9vMmp7F8wx4+LCxn4ekjrC5HqbDmcDgYN26c1WWEvLBsxpk5ehhJcZGsc5ZaXYpSSg2KsAx7u02YNzmNjYXltLjcVpejlFIDLizDHuCi7DSq65vJ23fU6lKUUmrAhW3YfzszlUi7jXUF2pSjlAp9YRv2cVERnD0hmXXOUr3lSykV8sI27AEuzEpjb0Ude8prrS5FKaUGVFjeetlqflY6P/3rDtY7S5mYFm91OafO1QIt9dDSCM3eYVfTribPuu4WcDd7ht1NGxcYN7hdYIxnvHVe23z3iWUYH8Y5Md2qbdycPN5pOT4s67D8pNm+/i9O/7enLJAyCS791wHZdViH/YjEGLJPG8I6Zym3zZ1gdTknNDfAsQNQvR9qyqChGhqqvMMeXu4W/9YhNhA72OwnxsUGIu3mtc4XzzjiHRfvuK3rcehhnA7z2wrqUF8PyzotP2mBj8evXWCrQdY0cK0MYR324GnK+Y8NRRytbWJYXOTgvGntETi6zxPm1SUngr26xPOq7eZJQkcsRA898YpLheSJnvGoIZ7ljmiIiIaIKIiI8Qwd3mGEd5k9EuwOsEV4Xu3HW6fFDrawbuVTKqRo2Gen85sPithQWMZ3zuj7I8i9ammEQ19AyTbP60AeVH1z8jqOOEjMgKGj4LTTPcOh3un44RCT6AnziEH6MlJKhRyfwl5ELgH+HbADvzfGPNlh+e3AHYALqAFuNcYUeJc9CPzAu+wuY8y7/iv/1OWMGEpaQhTrnKWnHvbGQNU+KMnzvrbB4S88beTgCfCR34Izb/WckQ8d5XlFJ2qTgVJqQPUa9iJiB5YDFwElwDYRWdMa5l4vG2Oe866/CPgVcImIZAPXAlOBEcA6EZlkjHH5+Tj6zWYT5mel8+bnB2lscREVYe/bDozxBPunL8Kud080wUTEwMgz4Kz/A6NmwchcGHKa3+tXSilf+HJmfyZQZIwpBhCRV4HFQFvYG2OOtVs/jhO3MiwGXjXGNAJfi0iRd3//8EPtfnNhVhqvbP2GT4orOW9Sqm8b1R+FL1ZC/h+hbIenKWbK5TB6tifc06aCPexbyZRSAcKXNBoJ7G83XQLM7riSiNwB3ANEAhe02/bjDtuO7FelA2jOxBSiHTbWO0t7Dntj4JuPIf9FKFgNLQ0wYiZc8QxMuwqiEgatZqWU6gtfwr6rxuRONyEbY5YDy0XkeuBhYImv24rIrcCtAKNHj/ahJP+Kdtg5d2Iq65xlPLbIIB3bz+sq4fNXPGfxRwohMgFm3ADfWuK5oKqUUgHOl7AvATLaTY8CDvaw/qvAs33Z1hjzPPA8QG5uriVPs1yUncY6Zyk7Dx8n67QhnpnHDsK6n8GOv3guso6aBYuXw9R/gsg4K8pUSql+8SXstwGZIjIOOIDnguv17VcQkUxjzG7v5OVA6/ga4GUR+RWeC7SZwFZ/FO5v86akAbCuoJSs4Qmw/SX42794Qv5bN3vO4tOnWlylUkr1T69hb4xpEZE7gXfx3Hr5gjFmh4gsA/KMMWuAO0XkQqAZOIqnCQfveivxXMxtAe4IpDtx2ktLiGZGRiLbd+yAQ0uhaB2MPgcW/wckB9DTtUop1Q8SaD0+5ubmmry8vMF/Y2P44OWnmbXrV8Q7QC5aBrP+tz5FqpQKCiKSb4zJ7W653hsInq4L3ryLC4o38nd3NhXn/hsLZ8+xuiqllPKb8A57txvyX4D3HwXAXP4r7l83hin7olhocWlKKeVP4Rv2lV/Dmh/B3s0wfh4s+g2SOJr5B7/itbz9NDS7iHb08WlapZQKUOHXIG0MbPtvePYcOPQ5LPwN3PQGJHru778wO52GZjdbio5YXKhSSvlPeIV9/VFY+T1Yew+MPgt++A/PLZXtHqKaPS6Z+KgI1jn1t2mVUqEjfJpx9m+FVT+A4wfhop/D2Xd2eadNZISNuZNSWe8sw+022GzaG6VSKviF/pm92w2bfwkvXOI5g7/lPZhzV4+3VM7PSqPseCNfHawexEKVUmrghPaZ/fFSeONWKN4IU78DC5/x/KpTL+ZNTsMmnqdpp49KHPg6lVJqgIXumX3RenhuDnzzieci7FUv+BT0AMPiIskdk8Q6Z9kAF6mUUoMj9MLe1ey5b37FdyA2BW7d0OkirC8uzE6j4NAxDlTVD1ChSik1eEIr7I/u9bTNb3nG03nZrRsgLatfu5qflQ7AB3pXjlIqBIRO2B8pgufOgyO74eoXPe3zjph+725CajzjUuJ4X5tylFIhIHTCPmk85N4Mt2/y9DfvBxdmpfHxngpqGlv8sj+llLJK6IS9zQYX/QyGjfXbLudnpdPkcrNem3KUUkEudMJ+AOSOGcb41Dj+ff1uml1uq8tRSql+07DvQYTdxr9cmkVxeS2vbP3G6nKUUqrfNOx7MT8rjbPHJ/Pr93dRXd9sdTlKKdUvGva9EBEeujyLqvpm/nNDkdXlKKVUv2jY+yBn5FC+e8Yo/rBlL/sr66wuRyml+kzD3kf3LpiM3SY8+bedVpeilFJ9pmHvo+FDo7n1vPGs/eIQ+fsqrS5HKaX6RMO+D26bO560hCh+sdaJMcbqcpRSymca9n0QGxnBvRdP5rNvqnjri0NWl6OUUj7TsO+j754xiqzThvDkOztpaHZZXY5SSvlEw76P7Dbh4cuzOFBVz4t/32t1OUop5RMN+36YMzGF+VPSWP5BERU1jVaXo5RSvdKw76cHL8uirtnFM+t2W12KUkr1SsO+nyamxXPD7NG8vPUbisqOW12OUuXqZecAAAn6SURBVEr1SMP+FNw9P5NYh53/97Y+aKWUCmwa9qcgOT6KOy+YyAc7y/ho9xGry1FKqW5p2J+iJeeMZdSwGH6xtgCXWx+0UkoFJg37UxTtsLP00insPHycVfn7rS5HKaW6pGHvB5dPO40zRifyxDs7+epAtdXlKKVUJxr2fiAiPPO/ZhIXGcF1z3+sHaUppQKOT2EvIpeISKGIFInI0i6W3yMiBSLyhYisF5Ex7Za5RGS797XGn8UHktHJsfz59rNJSYjixt9vZUuRXrBVSgWOXsNeROzAcuBSIBu4TkSyO6z2GZBrjJkOrAKeares3hgzw/ta5Ke6A9KIxBheu+0sRifFcvOL21hXUGp1SUopBfh2Zn8mUGSMKTbGNAGvAovbr2CM2WCMaf0Jp4+BUf4tM3ikJUTz6q1nMWV4ArevyOfNzw9aXZJSSvkU9iOB9reZlHjndecHwDvtpqNFJE9EPhaRK7vaQERu9a6TV15e7kNJgW1YXCQv/e/ZnDF6GHe/+hkr8/QuHaWUtXwJe+liXpc3lIvIjUAu8HS72aONMbnA9cAzIjKh086Med4Yk2uMyU1NTfWhpMCXEO3gj7ecyZyJKdy/6gte3PK11SUppcKYL2FfAmS0mx4FdGqbEJELgYeARcaYtq4gjTEHvcNiYCMw8xTqDSoxkXZ+vySXBdnpPPZmAcs3FFldklIqTPkS9tuATBEZJyKRwLXASXfViMhM4Hd4gr6s3fxhIhLlHU8B5gAF/io+GERF2Fl+wxksnjGCp98t5Ol3d+pPGiqlBl1EbysYY1pE5E7gXcAOvGCM2SEiy4A8Y8waPM028cCfRQTgG++dN1nA70TEjeeL5UljTFiFPYDDbuNX18wgNtLO8g17qG108cgV2dhsXbWQKaWU/0mgnWXm5uaavLw8q8sYEMYYfv6Wkxe2fM20kUP58YWZXDAlDe8XpFJK9ZuI5Huvj3ZJn6AdRCLCT6/I4pdXn05VfRM/+GMeVy7fwobCMm3aUUoNKD2zt0izy81fPi3hN+uLOFBVz8zRifzkwkl8OzNFz/SVUn3W25m9hr3FmlrcrMov4T8+2M3B6gZyxwzjJxdN4pwJyRr6SimfadgHicYWFyvzSlj+QRGHjzVw5rgkfnLhJM6ekGx1aUqpIKBhH2Qaml28uvUb/nPjHsqON3J6RiLzJqfy7cwUpo9KxGHXyyxKqc407INUQ7OLlz/5htXbD/DlgWqMgfioCM4an8y3M1M4NzOF8Slx2tSjlAI07ENCVV0Tf99TwebdR/ioqJz9lfUAjBgazZyJnuCfMzGFlPgoiytVSllFwz4EfVNRx+aicj7afYS/76mgur4ZgNSEKMYlxzEmOZaxKXGMbTceH9Xr83NKqSCmYR/iXG7DVweq+UdxBXvKathXUcfXFbWUH288ab2U+CjGpcQyJjmOUcNiGBYbSWKsg6S4yJPGYxx2bRpSKgj1FvZ6uhfk7Dbh9IxETs9IPGl+bWML+yrq2FtR63kdqWVvRR2bdpVT1uGLoL3ICBvDYh0Mi/V8CcRF2Yl22Ilx2ImJ9L7aT3uHURF2HHbBYbcRYRMcETYcNhuOCCHCZiPSbiPCLp6XzYZNwGYT7CLYbYKtbYh+2Sg1ADTsQ1RcVATZI4aQPWJIp2UtLjdV9c1U1TVxtK6ZytqmtvGjdU0crfWMV9U1cai6mfpmF/VNrrZhY4t7QGsXAbsINpsggE0EEe/Qu7zjMvAMpd0+pN289l8graNtQ+9WJ6Zbl3f+0un2a8jH76e+fI0Fy5decFQZHKacNoTfXjcwHQNr2IehCLuNlPiofl/QdbkNDc0nwr+h2UVDs5tmt5vmFjctbkOTy02Ly9DicreNN7vcNLsNbrfB5Ta4jWfoMq3zwG1Ono/nD2638QyNwRhPP0NuAwbv0EDrzywY77TBuy6tyz3zOLFq2w8ztDZnnpjufNzdNXj62hTapwbTwGpd7ZYJlkKDRMawmAHbt4a96jO7TYiLiiBOL/oqFTT0CR2llAoDGvZKKRUGNOyVUioMaNgrpVQY0LBXSqkwoGGvlFJhQMNeKaXCgIa9UkqFgYDrCE1EyoF9p7CLFOCIn8oJBKF2PBB6xxRqxwOhd0yhdjzQ+ZjGGGNSu1s54ML+VIlIXk89vwWbUDseCL1jCrXjgdA7plA7Huj7MWkzjlJKhQENe6WUCgOhGPbPW12An4Xa8UDoHVOoHQ+E3jGF2vFAH48p5NrslVJKdRaKZ/ZKKaU60LBXSqkwEDJhLyKXiEihiBSJyFKr6/EHEdkrIl+KyHYRCbpfYReRF0SkTES+ajcvSUTeF5Hd3uEwK2vsq26O6TEROeD9nLaLyGVW1tgXIpIhIhtExCkiO0Tkbu/8oPycejieYP6MokVkq4h87j2mn3nnjxORT7yf0WsiEtnjfkKhzV5E7MAu4CKgBNgGXGeMKbC0sFMkInuBXGNMUD4MIiLnATXAn4wxOd55TwGVxpgnvV/Kw4wxD1hZZ190c0yPATXGmH+zsrb+EJHTgNOMMZ+KSAKQD1wJfJ8g/Jx6OJ5rCN7PSIA4Y0yNiDiAj4C7gXuAvxhjXhWR54DPjTHPdrefUDmzPxMoMsYUG2OagFeBxRbXFPaMMZuAyg6zFwN/9I7/Ec8/xKDRzTEFLWPMIWPMp97x44ATGEmQfk49HE/QMh413kmH92WAC4BV3vm9fkahEvYjgf3tpksI8g/YywDviUi+iNxqdTF+km6MOQSef5hAmsX1+MudIvKFt5knKJo8OhKRscBM4BNC4HPqcDwQxJ+RiNhFZDtQBrwP7AGqjDEt3lV6zbxQCXvpYl7wt0/BHGPMGcClwB3eJgQVeJ4FJgAzgEPAL60tp+9EJB54HfixMeaY1fWcqi6OJ6g/I2OMyxgzAxiFpyUjq6vVetpHqIR9CZDRbnoUcNCiWvzGGHPQOywD3sDzIQe7Um+7amv7apnF9ZwyY0yp9x+jG/gvguxz8rYDvw68ZIz5i3d20H5OXR1PsH9GrYwxVcBG4CwgUUQivIt6zbxQCfttQKb36nQkcC2wxuKaTomIxHkvMCEiccAC4KuetwoKa4Al3vElwF8trMUvWkPR658Ios/Je/HvvwGnMeZX7RYF5efU3fEE+WeUKiKJ3vEY4EI81yI2AFd5V+v1MwqJu3EAvLdSPQPYgReMMY9bXNIpEZHxeM7mASKAl4PtmETkFeB8PF2xlgKPAquBlcBo4BvgamNM0Fzw7OaYzsfTPGCAvcBtre3dgU5EzgU2A18Cbu/sf8HTzh10n1MPx3MdwfsZTcdzAdaO5wR9pTFmmTcjXgWSgM+AG40xjd3uJ1TCXimlVPdCpRlHKaVUDzTslVIqDGjYK6VUGNCwV0qpMKBhr5RSYUDDXimlwoCGvVJKhYH/D4LfMAdoXwWTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#进行可视化\n",
    "x = np.arange(30)\n",
    "plt.plot(x,value1,label='cheerful')\n",
    "plt.plot(x,value2,label='so-so')\n",
    "plt.plot(x,value3,label='sad')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "z2OC3tyoJNeN"
   },
   "source": [
    "可以发现，从10轮左右开始，我们的状态概率分布就不变了，一直保持在  \n",
    "[0.23076934,0.30769244,0.4615386]\n",
    "\n",
    "### https://zhuanlan.zhihu.com/p/37121528"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3Taf1Z8nJU12"
   },
   "source": [
    "#### M-H采样python实现  \n",
    "https://zhuanlan.zhihu.com/p/37121528\n",
    "\n",
    "假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 $Q(i,j)$ 的条件转移概率是以 $i$ 为均值,方差1的正态分布在位置 $j$ 的值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 323
    },
    "colab_type": "code",
    "id": "eUoFG0tYJx9e",
    "outputId": "b404de1a-535b-4772-f90d-91f1f17ebe1e"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXiU5b3/8fc3YQmi7CCyGZRdgoFGFlEWxUqrCLVa0FK1Vam16rG2tLFyiXKwh1P6U9ujtqXUA62oKLUpVqsFoVgVKghRNlGwCAEFBEF2SLh/f8wkJ8skc0+S2T+v68rlzDP38+QbnPk+99yrOecQEZH0khHvAEREJPaU/EVE0pCSv4hIGlLyFxFJQ0r+IiJpqEG8A6isTZs2Ljs7O95hiIgklXfeeecz51xb3/IJl/yzs7NZtWpVvMMQEUkqZvZxJOXV7CMikoaU/EVE0pCSv4hIGkq4Nv9QTp48SVFREceOHYt3KJIEsrKy6NSpEw0bNox3KCIJKymSf1FREWeccQbZ2dmYWbzDkQTmnGPv3r0UFRXRtWvXeIcjkrCSotnn2LFjtG7dWolfwjIzWrdurW+JImEkRfIHlPjFm94rIuElTfIXEZH6kxRt/lWMGVO/13vxxbBFHnroIZ5++mkyMzPJyMjgt7/9LYMGDarfOMoZMWIEv/jFL8jLy6v1NebMmcPkyZPp1KkThw4d4pxzzmHq1KlceOGFANx///0MGzaMUaNGhTy/oKCAHj160KdPn5Cv/+Y3v+G0007jhhtuiDje/fv38/TTT3P77bcDsHPnTu666y4WLFhQi79URCKVnMk/xpYvX85f//pXVq9eTePGjfnss884ceJEvMPyMn78eB577DEAli5dytVXX83SpUvp3bs306ZNq/HcgoICrrzyypDJv7i4mNtuu63Wce3fv58nnniiLPl36NBBiV8khtTs4+GTTz6hTZs2NG7cGIA2bdrQoUMHAKZNm8YFF1xA3759mTRpEqU7o40YMYIf/OAHDBs2jN69e7Ny5UquvvpqunfvzpQpUwDYunUrvXr14sYbb6Rfv35cc801HDlypMrv//vf/86QIUMYMGAA1157LYcOHQIgPz+fPn360K9fP370ox+F/TtGjhzJpEmTmDVrFgA33XRTWcKtfK233nqLhQsXMnnyZHJzc9myZQsjRozgpz/9KcOHD+eXv/wlDzzwAL/4xS/Krv/UU09x4YUX0rdvX95++22AKmX69u3L1q1byc/PZ8uWLeTm5jJ58mS2bt1K3759gUAH/7e//W1ycnLo378/S5cuBQLfZK6++mpGjx5N9+7d+fGPf+z7v1BEKlHy9/DlL3+Z7du306NHD26//XaWLVtW9todd9zBypUrWbduHUePHuWvf/1r2WuNGjXi9ddf57bbbmPs2LE8/vjjrFu3jjlz5rB3714ANm3axKRJk3jvvfdo1qwZTzzxRIXf/dlnnzF9+nQWL17M6tWrycvL4+GHH2bfvn38+c9/Zv369bz33ntlN5RwBgwYwPvvv1/hWKhrXXjhhVx11VXMnDmTwsJCzj33XCBQY1+2bBk//OEPq1z78OHDvPXWWzzxxBN85zvfqTGOGTNmcO6551JYWMjMmTMrvPb4448DsHbtWp555hluvPHGstE7hYWFzJ8/n7Vr1zJ//ny2b9/u9XeLSEVK/h5OP/103nnnHWbNmkXbtm0ZP348c+bMAQJNKYMGDSInJ4clS5awfv36svOuuuoqAHJycjjvvPM466yzaNy4Meecc05Z0urcuTNDhw4FYOLEibzxxhsVfveKFSvYsGEDQ4cOJTc3l7lz5/Lxxx/TrFkzsrKyuOWWW3jhhRc47bTTvP6WUHs2R3Kt8ePHV/vaddddB8CwYcP44osv2L9/v1dMlb3xxht861vfAqBXr16cffbZfPDBBwBceumlNG/enKysLPr06cPHH0e0lpWIBKnN31NmZiYjRoxgxIgR5OTkMHfuXCZMmMDtt9/OqlWr6Ny5Mw888ECF8eWlzUQZGRllj0ufFxcXA1WHJVZ+7pzjsssu45lnnqkS09tvv81rr73Gs88+y2OPPcaSJUvC/h1r1qyhd+/eFY41aNDA+1pNmzat9tqh/pYGDRpw6tSpsmM+4+9D3aBKlf93zMzMLPt3FJHIqObvYdOmTXz44YdlzwsLCzn77LPLElmbNm04dOhQrTost23bxvLlywF45plnuOiiiyq8PnjwYN588002b94MwJEjR/jggw84dOgQBw4c4Ktf/SqPPvoohYWFYX/XsmXLmDVrFrfeemuF49Vd64wzzuDgwYPef8v8+fOBQM29efPmNG/enOzsbFavXg3A6tWr+fe//x322sOGDWPevHkAfPDBB2zbto2ePXt6xxEXY8aE/hFJUMlZ8/cYmlmfDh06xJ133sn+/ftp0KAB3bp1Y9asWbRo0YJbb72VnJwcsrOzueCCCyK+du/evZk7dy7f/e536d69O9/73vcqvN62bVvmzJnDddddx/HjxwGYPn06Z5xxBmPHjuXYsWM453jkkUdCXn/+/Pm88cYbHDlyhK5du/KnP/2pSs3/4MGDIa81YcIEbr31Vn71q1953dhatmzJhRdeyBdffMGTTz4JwNe//nX+8Ic/kJubywUXXECPHj0AaN26NUOHDqVv37585Stf4fvf/37ZdW6//XZuu+02cnJyaNCgAXPmzKlQ4xeRurOavmLHQ15enqu8mcvGjRurJKxUsHXrVq688krWrVsX71BSTszfM9XV8mNcUZH0ZWbvOOe8JwZ51fzNbDTwSyATmO2cm1Hp9XuAW4BiYA/wHefcx8HXbgRKh6JMd87N9Q1OJJFk579U7WuzN+4CYFTvM2MVjkidhE3+ZpYJPA5cBhQBK81soXNuQ7lia4A859wRM/se8HNgvJm1AqYCeYAD3gme+3l9/yHJKDs7W7X+BFdTwg9lcfAmUCr03GmR+POp+Q8ENjvnPgIws2eBsUBZ8nfOLS1XfgUwMfj4cmCRc25f8NxFwGig6tAVkQQSadIPd52tM66ol+uJ1Bef0T4dgfIzaYqCx6pzM/C3SM41s0lmtsrMVu3Zs8cjJJHoKFizo94Sf3nRuKZIXfgk/1Dr44bsJTaziQSaeEqnbHqd65yb5ZzLc87ltW3b1iMkkfqXnf8Sd88PP2S2LtcvWLMjatcXiYRP8i8COpd73gnYWbmQmY0C7gOucs4dj+RckXiLVc387vmF+hYgCcEn+a8EuptZVzNrBEwAFpYvYGb9gd8SSPy7y730KvBlM2tpZi2BLwePJZW9e/eSm5tLbm4u7du3p2PHjmXPo7W65+rVq3nllVdCvrZ48WKaN29O//796dGjB8OHD+fll18ue/3xxx8vmyQVypIlS1ixYkW1r//5z38uW29n4sSJFBQUeMd96tQpZsz4v8FgJSUlXHzxxd7nx0M8krFuABJvYTt8nXPFZnYHgaSdCTzpnFtvZtOAVc65hQSaeU4Hng9O8d/mnLvKObfPzP6TwA0EYFpp528yad26ddms1wceeIDTTz/daxXNUiUlJWRmZkb0O1evXs26desYPXp0yNdHjhxZlpRXr17N1772Nf7whz8wfPjwChOmQlmyZAlt2rRh8ODBVV4rLi7ma1/7WkSxllea/PPz84HAEgz//Oc/a329aItnEs7Of0kdwRI3Xss7OOdeds71cM6d65x7KHjs/mDixzk3yjl3pnMuN/hzVblzn3TOdQv+/G90/oyKCtbsYOiMJXTNf4mhM5ZEtZ11zJgxfOlLX+K8885j9uzZQCCBtmjRgilTpjBw4EDefvttFi5cSM+ePbn44ou58847GTduHBCYPXzTTTcxcOBA+vfvz4svvsjRo0eZNm0a8+bNIzc3N+zs2gEDBnDfffeVrds/ZcoUHn30UQAeeeQR+vTpw/nnn8/EiRPZsmULs2fPZubMmeTm5vLWW28xceJEfvjDHzJy5Eh++tOfMnv2bO6+++6y67/66qtcfPHF9OjRg7/9LdCXX7nM6NGjeeONN8jPz+fgwYPk5uZyww03lP1bQODGcM8999C3b19ycnLK/q7Fixdz6aWXcvXVV9OzZ09uuOGG+vhfE1ZtEn8DC4zcqfxT2/H9+gYg8ZKcyzvUoGDNDu59YS1HT5YAsGP/Ue59YS0A4/rXNEipdubOnUurVq04cuQIeXl5fP3rX+eMM87gwIEDDBgwgOnTp3PkyBF69OjBm2++SZcuXfjGN75Rdv60adMYPXo0c+bM4fPPP2fQoEG899573H///axbt64siYczYMAA/ud//qfK8Z///Od8/PHHNGrUiP3799OiRQtuueUW2rRpU5a8n3jiCbZs2cJrr71GRkZG2U2s1Pbt21m2bBkffvgho0aNKltnKJQZM2Ywe/bssm9K5Rdee/7559mwYQPvvvsue/bs4YILLmDYsGFA4NvLhg0baNeuHYMHD2bFihUhv5nUl9ok3XC19NIbQOWx/j6xlF1bM4UlRlJuYbeZr24qS/yljp4sYearm6Ly+x555BHOP/98hgwZQlFREVu2bAECa/mXNp9s2LCBnj17cvbZZ2NmZUsfQ2Cjloceeojc3FxGjhzJsWPH2LZtW8RxVLdMx3nnncfEiROZN28eDRs2rPb8a6+9loyM0G+Hb3zjG2RkZNCzZ086d+5cYZG7SLzxxhtcf/31ZGZm0r59ey666CJKl/IYPHgwZ511FpmZmeTm5rJ169Za/Q4fkSb+0tq9r9p8C9AoIIm1lKv579x/NKLjdbF48WJef/11VqxYQZMmTbjooovKVvps0qRJ2RLHNa2f5JyjoKCgbLOUUq+//npEsYRaqhkCTTbLli3jL3/5C9OnT692RnG6LNVcOfHPXvBgyHK3XDMVqP3krEhvAHfPL4zKN1OR6qRczb9DiyYRHa+LAwcO0KpVK5o0acL69etZuXJlyHLnnXcemzZtYvv27TjnypY+Brj88sv51a9+VfZ8zZo1QGTLKRcWFvKzn/2sSkdvSUkJRUVFXHLJJcycOZM9e/Zw5MiRiJdqfv7553HO8cEHH7B9+3a6d+9OdnY2a9aswTnH1q1beeedd4DA3gBAyOQ9bNgwnn32WUpKSti1axdvvvlmnTaoj1Sv+14OX6ic+uiMjeQaav+XWEq55D/58p40aVhxZE2ThplMvrz+14O/4oorOHLkCOeffz7Tpk1j0KBBIcuddtppPPbYY4waNYqLL76YDh060Lx5cwCmTp3KkSNHynb7euCBBwC45JJLePfdd+nfv3/IDt+lS5fSv39/evbsyV133cUTTzzB8OHDK5QpLi7m+uuvp1+/fgwYMICf/OQnZUtBP/fcc/Tv35+33nor7N/ZrVs3hg0bxpgxY5g1axaNGjVi+PDhdOzYkZycHPLz88nNzS0rf/PNN9OvX78qHbfXXHMNvXr14vzzz2fUqFE8/PDDtGvXLuzvry/HSvxXsK3PUTiRXCvS/gKR2krJJZ0L1uxg5qub2Ln/KB1aNGHy5T3j/pX60KFDnH766Tjn+O53v0tOTg533nlnXGNKZZXfM9XVqkM1+4zqfWbkHaweHbU+NfvZCx7EgEsrNxupw1fCiMqSzslmXP+OcU/2lf36179m3rx5HD9+nLy8vCq7aUn0RNKcEs0lmbfOuKIslur6GqCatVNE6lnKNfskqsmTJ1NYWMjGjRv54x//SFZWVrxDSguDHlrkXTYWa/EPPbeVVzk1/0i0JU3N3zlXZdSJSCjlmzJ3HQwsv1FTTRtitwnLvFuH0NXzm8hrG3dVbf4RqSdJUfPPyspi7969NQ4VFIFA4t+7dy9ZWVnezT19OzSPclQV/duzA1jvdommpKj5d+rUiaKiIrTWv/jIysrioX986l2+ffMoNsFV0xE8qveZXk07izfu0taQEhVJkfwbNmxI165d4x2GJJHFv90QvhDx3XO3YYZx8lT4+v2Sjbu4JAbxSHpJiuQvUkF1wyoBXnyRbvf6NffEu0Y9vGc7r9r/qbAlRCKn5C8pp9ijsbxp4zBv/RgtsObb/KPln6W+KflLSlncbSCzwxdjyDmtox6Lr1ZNG7HvcHQ2BRKpTlKM9hHxsWzT7vCFiH9zT2UDurT0Kqe1f6Q+KflLyvDpPE1UvjekSBenE6mOmn0kJfjOiK1zrb+mzuYYiGRxOpGaqOYvaSPWk7ki5XtjUvOP1Aclf0l6vrX+qE7mqieJ1h8hqUvJX9JCqiVV1f6lrpT8Jaml4uqXqXajksSk5C8pLxmTadhJaKj2L3Wj5C9JKxVr/aUSaRKapCYN9ZSUloy1/lCq248gm/rdb1jSh2r+kpR8av2dWjaJQSTRkyo3LklMqvlLyurVvlm8Q6izxg0yOF5c/bqesxc8yOIFD1a9UWjDdwlDNX9JOj61/lZNG8Ugkui7uHvbeIcgKUrJX1KS72JpycCn+SqVO78lOpT8Jan4DG9MlVp/qVRovpLEo+QvKSeVav2lfNYlUu1fIqHkL0nDp9afqiNkkmFdIkkuSv4iSUK1f6lPSv6SFNK51l+qffMsfWCl3ni9l8xstJltMrPNZpYf4vVhZrbazIrN7JpKr5WYWWHwZ2F9BS6Sji7xuMGp9i8+wk7yMrNM4HHgMqAIWGlmC51zG8oV2wbcBPwoxCWOOudy6yFWSVM+tf5E36hFJNH41PwHApudcx85504AzwJjyxdwzm11zr0HVD8VUSSK0qlD1Gfcv1b8lHB8kn9HYHu550XBY76yzGyVma0ws3ERRSdpb9BDi8KWSbdav8b9S33wSf4W4lgku0h3cc7lAdcDj5rZuVV+gdmk4A1i1Z49eyK4tKS6XQdPhC2TTrX+Uj7r/XdV7V9q4JP8i4DO5Z53Anb6/gLn3M7gfz8C/gH0D1FmlnMuzzmX17at1jKRgG73pt9sXl8+6/1HUkOT9OOT/FcC3c2sq5k1AiYAXqN2zKylmTUOPm4DDAU21HyWSECxR/ZKxdm8vnw+vN/83fKoxyHJKez7xzlXDNwBvApsBJ5zzq03s2lmdhWAmV1gZkXAtcBvzWx98PTewCozexdYCsyoNEpIJKSCNTvClsm0UC2S6cNn2OebW/bFIBJJRl7r+TvnXgZernTs/nKPVxJoDqp83ltATh1jlDR09/zCsGVG9moXg0gSW6YZJa7mr0hTCtYyfZw+hlKRJgxKwvGp9euNG+BzA3xqxbYYRCLJRp8hSTg+tX6fJg8RqZ6SvyQdvWkr8lnTyGe+hKQXfY4kofgkKdX6I+czX0LSi5K/JBQlqdrxqf37zJuQ9KHkLwmj130vhy2zdcYVMYgkNfnMm5D0oeQvCeNYibJTXfiscdRv6isxiESSgZK/JASftv5Hx2tl8Jr4rHH0xfGSGEQiyUDJXxKCT1v/uP6RLCabnhpmhJ/17DOPQlKfkr/EnU8yatY4MwaRJL/hPcNP+vKZRyGpT8lf4s4nGb334OgYRJIa9KEWH15r+4hE0+wFD4Y8fss1UwFokN7rt0XMZx5Edv5LGjmV5pT8Ja6y819idjWvld4URvU+E8b8JnZBiaQBJX+RFLR1/W9YvHFXleOl36YgMOxTzWnpS82DEjc+k7p8Zq5K7WjYZ3pT8pe40aSu6PKZ9DWlYG0MIpFEpOQvceEzqcsneUn1fCZ9aa3/9KXkL3HhM6nLJ3lJzTq1bBK2jJZ7Tk9K/hJzPk0NPklLwuvVvlnYMlpJNT0p+UvM+TQ1+CQtEak9JX9JOJrTVb98Rkxl52ut/3Sj5C8x5ZNkLtXwznqXabqlSkWa5CWSBkb2asfijbuqXUoDoGB8rlZOTSNK/hI7Y8YwO8Ss0/I0qSt+7p5fqOSfRtTsIzETarkBiR2fEVTf/N3yGEQiiUDJX2LCZyy5hndGl88Iqje37ItBJJIIlPwlJnzGkmt4Z/T5dPxqp6/0oOQvCaFxA70VY2FkL+30JQH6xEnU+czovbh72xhEIqAbrQToXSBRF25Gr8agx5bPjVaTvlKfkr9EVbd7wycRn6YIEalfSv4SVcVhluxvmKFafzz4zKfQsM/UpuQvUeMzamR4T9X6E5WGfaY2JX+JGo0aSWw+8yo07DN1KflLVPg0GWjUSXz5zKvQDTx1ea3tY2ajgV8CmcBs59yMSq8PAx4F+gETnHMLyr12IzAl+HS6c25ufQQuCWDMmNDHX3zRq8lAwzvjr3GDDI4XnwKoftG3GVfEMCKJlbBVLzPLBB4HvgL0Aa4zsz6Vim0DbgKernRuK2AqMAgYCEw1s5Z1D1uSnWr9iUHDPtOXT81/ILDZOfcRgJk9C4wFNpQWcM5tDb52qtK5lwOLnHP7gq8vAkYDz9Q5cklYi7sNZHaYMqr1Jw4DwgzKkhTkU/3qCGwv97woeMyH17lmNsnMVpnZqj179nheWkTqg8/mOar9px6f5B9qILZvRcHrXOfcLOdcnnMur21b1QiT2T8/DH/z1pr9IvHnk/yLgM7lnncCdnpevy7nShIq7TysjqZ0JSafG/JlD/8j+oFIzPgk/5VAdzPramaNgAnAQs/rvwp82cxaBjt6vxw8Jilo2abdYctof97k9eHuw/EOQepR2OTvnCsG7iCQtDcCzznn1pvZNDO7CsDMLjCzIuBa4Ldmtj547j7gPwncQFYC00o7fyX1nDylbsNkpua49OI1zt859zLwcqVj95d7vJJAk06oc58EnqxDjJIEfLZo7NuheQwikWjKzn+JrRr3nxI02Fpipn3zrHiHIGG0atoobBkt+JYalPylzt7/9IuwZbQ/b3IY0CX8HEwt+JYalPylzoo+Pxq2jPbnTR4+s6+14FvyU/KXOvGp9Tdt7NW1JAnCZ/a1FnxLfkr+Uic+tf4h57SOQSRSn3w22fHZm1kSl5K/1JpPrV9vsOTks8lOuL2ZJbHpsym15lPrv0Rjx5OWav+pTclfokZt/clNtf/UpuQvtdLVY5VHtfUnP59x/5KclPylVsIt5KDNWlKDz7h/LfecnPQJlYj5fNi1WUvqyMrUWqypSMlf6l2mKVmkkvcf+mrYMur4TT5K/hIRn1r/yF7hOwoluYS7navjN/ko+YtIWP/2WMmz39RXYhCJ1Bclf/E26KFFYctoTfjU1SBM9f+L4yWxCUTqhZK/eNt18ESNr+vNlNo2/1f42r/PEGBJDJqFI158av2azZuixowpezi73KY9t1wztUpR7eWWPFRZEy+q9Qtoo/dUos+shKVav0RCG70nByV/qVHBmh1ha/2azZtefPZiVu0/8anNX2rks2mHZvOml/bNs9h/tOYKgWr/iU9VNqmWz0bdegOlJ59tObXmT2LTZ1eq5bNRt9r609fQc1uFLaO9fhOXkr+E5NNmq/X609u8W4eELaO9fhOXPr3yf8qN5/5JmPHcoPX6Bbq3axq2ff+bv1vudaOQ2FLylyqWf7S3wvPZCx6sUqZxgwxQR2/aW3TPiLBt+z7NhxJ7avaRKg4fLw5bRiN8pFSzxplhy2jRt8Sj5C8VVK71h+Kzsbekj/ceHB22jBZ9SzxK/lKBT63fZ2NvSS+Pjs8NW0ZDPxOLkr+UWVyuk7c6WrJZQhnXv6NXOQ39TBxK/gL4De3U9oxSk4mDu4Qto6GfiUPJXwC/6fjanlFqMn1cTtjtHkHr/iQKJX/x2oCjVdNGMYhEkp3Pdo9a9ycxaJx/mptSsNZrA44BXVpGPRZJMuUmBZY38eafhd3Qvd/UV7xGCUn0qOaf5sJ9SEFvEonM9HE5Ycto6Gf8eX2uzWy0mW0ys81mlh/i9cZmNj/4+r/MLDt4PNvMjppZYfDnN/UbvtSFb9urFm+TSPks+qahn/EVttnHzDKBx4HLgCJgpZktdM5tKFfsZuBz51w3M5sA/DcwPvjaFudc+EHAEjvBr+s/8RjaqbZ+qY15tw5h0EOLwm4EdNnD/2DRPSNiE5RU4FPzHwhsds595Jw7ATwLjK1UZiwwN/h4AXCpmcYFJrIlHonfUFu/1N6/7rssbJkPdx/W2P848enw7QhsL/e8CBhUXRnnXLGZHQBKl3zsamZrgC+AKc65f1b+BWY2CZgE0KVL+LHCUjefHjjGKY9yl6q5R2qjXEfw3E27OXkqMKSgutVh755f6D1JTOqPT80/VA2+8gCR6sp8AnRxzvUH7gGeNrMqWwA552Y55/Kcc3lt22rBsGhbt/NA2DKaySv1wXcpkEEPLYpyJFKZT/IvAjqXe94J2FldGTNrADQH9jnnjjvn9gI4594BtgA96hq01J5PJ5tm8kp96tSySdgyuw6e8No2VOqPT7PPSqC7mXUFdgATgOsrlVkI3AgsB64BljjnnJm1JXATKDGzc4DuwEf1Fr3UrNI47KXv72a2Cz+qXzN5pT71at+MPQePhy2ndf9jK2zN3zlXDNwBvApsBJ5zzq03s2lmdlWw2O+B1ma2mUDzTulw0GHAe2b2LoGO4Nucc/o/HCclHolfdX6JBt/9H7T0Q+x4zfB1zr0MvFzp2P3lHh8Drg1x3p+AP9UxRqkHr3mM7gF18kr0bJ1xRdhmRy39EDuavJkGlr6/22sJB3XySrRp16/EobV9Utzyj/Z6Nff4dMqJ1MmYMbxH1X0jKg8B/eJ4CV3zX/JaJE5qTzX/FPbpgWNeO3NBoFNOJBZ8Zo07NPwz2pT8U5jPeH5Qc4/E1oAuLb2GE+86eEKzf6NIyT9F+Y6Z7tuheZQjEalqZK92NG4QPv1o56/oUfJPQb3ue9lrzHSnlk1o3zwrBhGJVOU7/FMdwNGh5J9iet33MsdKwnfwZpqpnV/izuebZ2kHsNQvJf8U45P4QbN4JTG0b57lNbHQ4bfdqPhT8k8hvptjNG2sEb6SOHyHdDo0A7g+KfmnCN/E3zDDGHJO6/AFRWJoq+cN4MPdh7UAXD1R8k8BkWyH57vErkis+Wz9CIEF4DQEtO70/T+ZVFqlE4I7clWzSUZlGs8vCSn4vp5Hxdm/1W3+MnvBg7AA3m/ZpOKghRdfjGaUKUfJP4m9tnEXjuCHoQYNM0w1fkkKo3qfWfa+Dqfo86McOVGirUZrSc0+SWqx5wcE1NQjyeXS3md6LQEBsO/wCT49cCzKEaUmJUIMlEsAAAsUSURBVP8kVHlhrJpowTZJRgO6tKR7u6ZeZdftPMDqbZ9HOaLUo2afJPLpgWPe6/UANG6QoYlckrQWLf1/LCu3AXxN9h0+waCHFvGv+y6LQWSpQTX/JNFv6isRJf5WTRt5T58XSVTDe7bznpey6+AJzQOIgJJ/EsjOf4kvjpd4l+/bobk6wSRlDDmntXfz5Ye7D+sG4EnJP8FFOqW9VdNGWqxNUk4kzZcf7j5Mdv5LmgsQhpJ/Arvs4X94j+iBQOeuavySqiKdp3L3/EJtCFMDJf8ENKVgLdn5L3lvZp1B4IOhzl1JdZHeAHYdPEG3e7UgXChK/gmma/5LPLVim3d5Ay7RzF1JI6N6nxlR4ip2WhAuFCX/BHHZw/8gO/+liJp5jMCEGJF0c0nvM2mY4bMYdMCHuw9zjvoBKtA4/zgrWLOjylZ14ZZrgMAYfg3llHQ2vGe7iOa+nCLQD3D3/EImDu7C9HE50Q0wwSn5x9FlD//Du12/vL4dmmtEjwiBzWDaN89i6fu7IzrvqRXbeGrFNoae24p5tw6JUnSJTck/xr75u+Ve++tWR4lfpKqRvdrxx4L/5HjxqSqvVbc6KASWh87OfyktvwmYc5G0MkdfXl6eW7VqVbzDqHcFa3Zwz/xCqr41qwrV7KNmHpHw3v/0C4o+P+pVNtRN4dHxuYzr37H6k0Isqw4kxHLSZvaOcy7Pt7xq/lFW15o+BMbvaxinSHi92jejV/tmLP9oL4ePF0d8fmmfQMcWTZh8ec+abwRJTsk/SgrW7GDy84X8+rkH+XYtr9GqaSNN2hKphSHntK71DQBgx/6jZTeCVO0XUPKvJ6XJ/qRPu04YTRs30D67InVU+hn69MAx1u88ENEw6vJK+wUAZm/clTKfTyX/WipYs4MHX1zP50dO1ts1U+VNJZJISkcErd72OfsOn6jz9Q4fL66wp0anlk3oVeerxp6Sv4cpBWuZt2JbWc3BoNa1iFAyzOhzVjON4hGJotIm1NKbQHXzaWoaHRRK0edH+UpwgmYy9RUo+VdSsGYHM1/dxM79R+nQognZrZtU6LCtrzdMKbXri8RW6eetuo1ifCZZVlZ6lR37j/KD+YWs+ngfeWe3qpBLEu2mkLbJP1SSX/7RPsq/F3bsP8qO/X7DxiJ5wxjQUSN4ROKq8t7Wnx44xoZPvuBUHYe/OwKTyOav3M7JksC1duw/yr0vrC0rU/mmsOrjfTzzr+2UOEemGdcN6hz1eQdeyd/MRgO/BDKB2c65GZVebwz8AfgSsBcY75zbGnztXuBmoAS4yzn3ar1FX07lZD758p5A1X/kcf07UrBmB/e+sJajJwMbpESS5Gsjq2Em3dqermYdkQRW2jdQl1FC5ZUm/lJHT5bw4IvrOXbyVIXcU3n+T4lzZYs7RvMGEDb5m1km8DhwGVAErDSzhc65DeWK3Qx87pzrZmYTgP8GxptZH2ACcB7QAVhsZj2cc/7bUnkIlcwnL3gXHGVf68rfeWe+uqmsbH0345SnJh2R5DPknNZ8euAYm/cc4vjJknrt3ws1QKS6AYLP/Gt7fJM/MBDY7Jz7CMDMngXGAuWT/1jggeDjBcBjZmbB4886544D/zazzcHrLa+f8AN8k/nRkyVl3wTCqU27XylNyhJJbqXfAkqF20g+0kqkT/mSKK++4JP8OwLbyz0vAgZVV8Y5V2xmB4DWweMrKp1b7z0ePsm8fNkOLZp4N/NkZhjNmzTk88MnKoz2gUDbnmF0bJmlZC+Swir3DwBl3w6OnazakJFhgSxRPslnZhgZZpws8ZsMlGn+S1bXhk/yDxVB5VtSdWV8zsXMJgGTALp06eIRUkXlk3m45prStv/SZqKayjdpmMF/Xd2PkQnUQy8iiaF98Af8+xyBCk3UENhUJdTt4LpBnaMav0/yLwLKR9EJ2FlNmSIzawA0B/Z5notzbhYwCwILu/kGX6p8Mi/VMNMqtPkDNGmYWWG4VeXRPis++jymve0ikhrG9e8YchhndUM7E2G0T9hVPYPJ/APgUmAHsBK43jm3vlyZ7wM5zrnbgh2+VzvnvmFm5wFPE2jn7wC8BnSvqcO3tqt6RjLaR0Qk1dT7qp7BNvw7gFcJDPV80jm33symAauccwuB3wN/DHbo7iMwwodguecIdA4XA9+v75E+pSK984qIpDOt5y8ikgIirflrA3cRkTSk5C8ikoaU/EVE0pCSv4hIGlLyFxFJQ0r+IiJpSMlfRCQNKfmLiKQhJX8RkTSk5C8ikoaU/EVE0pCSv4hIGlLyFxFJQ0r+IiJpSMlfRCQNKfmLiKShhNvMxcz2AB+HKdYG+CwG4dS3ZI0bkjd2xR17yRp7ssYNgdibOufa+p6QcMnfh5mtimTHmkSRrHFD8sauuGMvWWNP1rihdrGr2UdEJA0p+YuIpKFkTf6z4h1ALSVr3JC8sSvu2EvW2JM1bqhF7EnZ5i8iInWTrDV/ERGpAyV/EZE0lNTJ38x+ZGbOzNrEOxZfZjbTzN43s/fM7M9m1iLeMdXEzEab2SYz22xm+fGOx5eZdTazpWa20czWm9l/xDumSJhZppmtMbO/xjuWSJhZCzNbEHyPbzSzIfGOyYeZ/SD4PllnZs+YWVa8Y6qOmT1pZrvNbF25Y63MbJGZfRj8b8tw10na5G9mnYHLgG3xjiVCi4C+zrl+wAfAvXGOp1pmlgk8DnwF6ANcZ2Z94huVt2Lgh8653sBg4PtJFDvAfwAb4x1ELfwSeMU51ws4nyT4G8ysI3AXkOec6wtkAhPiG1WN5gCjKx3LB15zznUHXgs+r1HSJn/gEeDHQFL1WDvn/u6cKw4+XQF0imc8YQwENjvnPnLOnQCeBcbGOSYvzrlPnHOrg48PEkhCHeMblR8z6wRcAcyOdyyRMLNmwDDg9wDOuRPOuf3xjcpbA6CJmTUATgN2xjmeajnnXgf2VTo8FpgbfDwXGBfuOkmZ/M3sKmCHc+7deMdSR98B/hbvIGrQEdhe7nkRSZJAyzOzbKA/8K/4RuLtUQIVm1PxDiRC5wB7gP8NNlnNNrOm8Q4qHOfcDuAXBFoRPgEOOOf+Ht+oInamc+4TCFR8gHbhTkjY5G9mi4Ptb5V/xgL3AffHO8bqhIm9tMx9BJom5sUv0rAsxLGk+qZlZqcDfwLuds59Ee94wjGzK4Hdzrl34h1LLTQABgC/ds71Bw7j0fwQb8H28bFAV6AD0NTMJsY3quhrEO8AquOcGxXquJnlEPif9K6ZQaDZZLWZDXTOfRrDEKtVXeylzOxG4ErgUpfYEy2KgM7lnncigb8OV2ZmDQkk/nnOuRfiHY+nocBVZvZVIAtoZmZPOeeSIRkVAUXOudJvWAtIguQPjAL+7ZzbA2BmLwAXAk/FNarI7DKzs5xzn5jZWcDucCckbM2/Os65tc65ds65bOdcNoE33IBESfzhmNlo4CfAVc65I/GOJ4yVQHcz62pmjQh0gi2Mc0xeLFAz+D2w0Tn3cLzj8eWcu9c51yn43p4ALEmSxE/wM7jdzHoGD10KbIhjSL62AYPN7LTg++ZSkqCjupKFwI3BxzcCfwl3QsLW/FPYY0BjYFHwm8sK59xt8Q0pNOdcsZndAbxKYATEk8659XEOy9dQ4FvAWjMrDB77qXPu5TjGlA7uBOYFKwsfAd+OczxhOef+ZWYLgNUEmmLXkMBLPZjZM8AIoI2ZFQFTgRnAc2Z2M4Gb2bVhr5PYrQ4iIhINSdfsIyIidafkLyKShpT8RUTSkJK/iEgaUvIXEUlDSv4iImlIyV9EJA39f6eBBVAiurkzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "\n",
    "def norm_dist_prob(theta):\n",
    "    y = norm.pdf(theta, loc=3, scale=2)\n",
    "    return y\n",
    "\n",
    "\n",
    "T = 5000\n",
    "pi = [0 for i in range(T)]\n",
    "sigma = 1\n",
    "t = 0\n",
    "while t < T - 1:\n",
    "    t = t + 1\n",
    "    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1,\n",
    "                       random_state=None)  #状态转移进行随机抽样\n",
    "    alpha = min(\n",
    "        1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))  #alpha值\n",
    "\n",
    "    u = random.uniform(0, 1)\n",
    "    if u < alpha:\n",
    "        pi[t] = pi_star[0]\n",
    "    else:\n",
    "        pi[t] = pi[t - 1]\n",
    "\n",
    "plt.scatter(pi, norm.pdf(pi, loc=3, scale=2), label='Target Distribution')\n",
    "num_bins = 50\n",
    "plt.hist(pi,\n",
    "         num_bins,\n",
    "         density=1,\n",
    "         facecolor='red',\n",
    "         alpha=0.7,\n",
    "         label='Samples Distribution')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RRUK6UCjKa4Y"
   },
   "source": [
    "#### 二维Gibbs采样实例python实现  \n",
    "\n",
    "假设我们要采样的是一个二维正态分布 $N(\\mu, \\Sigma)$ ，其中： $\\mu=(\\mu_{1}, \\mu_{2})= (5, -1)$ , $\\Sigma = \\begin{pmatrix}\n",
    "\\sigma^{2}_{1} &   \\rho \\sigma_{1}\\sigma_{2}b\\rho \\sigma_{2}& \n",
    "\\sigma^{2}_{2}\\end{pmatrix} = \\begin{pmatrix}\n",
    " 1& 1b1 & \n",
    "4\\end{pmatrix}$;\n",
    "\n",
    "而采样过程中的需要的状态转移条件分布为：\n",
    "\n",
    "$P(x_{1}|x_{2}) = N(\\mu_{1}+ \\rho \\sigma_{1}/\\sigma_{2}(x_{2} - \\mu_{2}), (1 - \\rho^{2})\\sigma^{2}_{1})$\n",
    "\n",
    "$P(x_{2}|x_{1}) = N(\\mu_{2}+ \\rho \\sigma_{2}/\\sigma_{1}(x_{1} - \\mu_{1}), (1 - \\rho^{2})\\sigma^{2}_{2})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 335
    },
    "colab_type": "code",
    "id": "3ZSOvaDcMHfv",
    "outputId": "0c3d6e18-ba4e-4e06-b61f-2585cc78715d"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEICAYAAAC+iFRkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAcOUlEQVR4nO3df5BV5Z3n8fcnrciOGn8AMYYGabWdCRpLMi0mZscQf2JmAjplNvijxJUUZSKjM6614mpIB8upaGbdSUoSZSaUmo1BYzbZ3hSGmCgzmx0x3cafoCwtorT4g4BKiID8+O4f57R7uNzbfbr79u3bfT6vqluc85znOfd7T1++97nPOfc5igjMzGxk+9BQB2BmZoPPyd7MrACc7M3MCsDJ3sysAJzszcwKwMnezKwAnOxtRJK0StK0oY7DrF442duwJGm9pLNLyq6Q9BuAiDgxIlb0so9JkkLSAYMYqlldcLI3GyT+ELF64mRvI1K25y9pqqQOSVslvSnpjrTav6b/viNpm6RPS/qQpJslvSLpLUn3SToss9/L022bJX2t5HlaJT0k6b9L2gpckT7345LekfS6pDsljcrsLyR9VdJaSX+QdIuk49I2WyU9mK1v1l9O9lYE3wa+HREfBo4DHkzLz0j/PTwiDomIx4Er0sfngGOBQ4A7ASRNBr4LXAocDRwGjC95rpnAQ8DhwA+BPcDfAWOBTwNnAV8taTMd+HPgU8B/BhanzzEBOAm4eACv3Qxwsrfh7Wdpj/kdSe+QJOJydgHHSxobEdsiYmUP+7wUuCMi1kXENuBGYFY6JHMR8L8i4jcR8T6wACidXOrxiPhZROyNiO0R8WRErIyI3RGxHrgb+GxJm9siYmtErAKeB36ZPv+7wMPAlPyHxKw8J3sbzi6IiMO7H+zfY+42BzgBeFFSu6S/6mGfHwNeyay/AhwAHJVu29C9ISLeAzaXtN+QXZF0gqSfS3ojHdr5e5JeftabmeXtZdYP6SFes1yc7G3Ei4i1EXEx8BHgNuAhSQezf68cYCNwTGZ9IrCbJAG/DjR2b5D074AxpU9Xsv494EWgOR1G+i+A+v9qzPrHyd5GPEmXSRoXEXuBd9LiPcAmYC/J2Hy3HwF/J6lJ0iEkPfEHImI3yVj8FySdnp40/Qa9J+5Dga3ANkl/Bnylai/MrA+c7K0IpgOrJG0jOVk7KyJ2pMMwtwL/Jx33/xSwBPgByZU6LwM7gL8BSMfU/wZYStLL/wPwFrCzh+e+HrgkrftPwAPVf3lmvZNvXmLWP2nP/x2SIZqXhzoes564Z2/WB5K+IOlP0jH/fwCeA9YPbVRmvXOyN+ubmSQncTcCzSRDQv56bHXPwzhmZgXgnr2ZWQHU3URNY8eOjUmTJg11GGZmw8qTTz75+4gYV2l73SX7SZMm0dHRMdRhmJkNK5Je6Wm7h3HMzArAyd7MrACc7M3MCqDuxuzNzIbKrl276OrqYseOHUMdSkWjR4+msbGRAw88sE/tnOzNzFJdXV0ceuihTJo0Can+JieNCDZv3kxXVxdNTU19authHDOz1I4dOxgzZkxdJnoASYwZM6Zf3zxyJXtJ0yWtkdQpaX4P9S5K76nZkim7MW23RtJ5fY7QzKyG6jXRd+tvfL0O40hqABYB5wBdQLuktohYXVLvUOAa4IlM2WRgFnAiyV1+fiXphIjY069ozcysX/KM2U8FOiNiHYCkpSSTQa0uqXcLcDvJ/N3dZgJLI2In8LKkznR/jw80cDOzwda6orW6+5tW3f31RZ5kP55976vZBZyWrSBpCjAhIn4u6fqStitL2o4vfQJJc4G5ABMnTswXuZkNC9mEOZTJrujyjNmXGyD6YKpMSR8C/hvwn/ra9oOCiMUR0RIRLePGVZzawcyGudYVrVXvLY807e3tnHzyyezYsYM//vGPnHjiiTz//PMD3m+enn0XMCGz3kgyl3e3Q4GTgBXpiYOPAm2SZuRoa2ZmGaeeeiozZszg5ptvZvv27Vx22WWcdNJJA95vnmTfDjRLagJeIznhekn3xoh4FxjbvS5pBXB9RHRI2g7cL+kOkhO0zcBvBxy1mdkItmDBAk499VRGjx7Nd77znarss9dkHxG7Jc0DlgMNwJKIWCVpIdAREW09tF0l6UGSk7m7gat9JY6ZWc+2bNnCtm3b2LVrFzt27ODggw8e8D5z/YI2IpYBy0rKFlSoO61k/Vbg1n7GZ2ZWOHPnzuWWW27h5Zdf5oYbbuDOO+8c8D49XYKZVd1IOQk7FFcP3XfffRxwwAFccskl7Nmzh9NPP51HH32UM888c0D7dbI3M6sjl19+OZdffjkADQ0NPPHEE720yMdz45iZFYCTvZlZATjZm5kVgJO9mVkBONmbmRWAr8Yxs6oYKZdbjlRO9mZmlbS21vf++sDDOGZmBeCevZnVnOe4r+xrX/saY8eO5dprrwXgpptu4qijjuKaa64Z0H7dszczqyNz5szh3nvvBWDv3r0sXbqUSy+9dMD7dc/ezKyOTJo0iTFjxvDUU0/x5ptvMmXKFMaMGTPg/TrZm5nVmS9/+cvcc889vPHGG1x55ZVV2aeHcczM6syFF17IL37xC9rb2znvvPOqsk/37M3MKhmiSyVHjRrF5z73OQ4//HAaGhqqss9cPXtJ0yWtkdQpaX6Z7VdJek7S05J+I2lyWj5J0va0/GlJd1UlajOzEWzv3r2sXLmSOXPmVG2fvSZ7SQ3AIuB8YDJwcXcyz7g/Ij4REacAtwN3ZLa9FBGnpI+rqhW4mdlItHr1ao4//njOOussmpubq7bfPMM4U4HOiFgHIGkpMJPkvrIARMTWTP2DgahahGZmBTJ58mTWrVtX9f3mGcYZD2zIrHelZfuQdLWkl0h69tmr/5skPSXpXyT9xYCiNTMbZBH13Vftb3x5kr3KPV+ZABZFxHHADcDNafHrwMSImAJcB9wv6cP7PYE0V1KHpI5Nmzblj97MrIpGjx7N5s2b6zbhRwSbN29m9OjRfW6bZxinC5iQWW8ENvZQfynwvTSwncDOdPnJtOd/AtCRbRARi4HFAC0tLfV5lM1sxGtsbKSrq4t67nSOHj2axsbGPrfLk+zbgWZJTcBrwCzgkmwFSc0RsTZd/UtgbVo+DtgSEXskHQs0A9UfjDKzITHSpjU+8MADaWpqGuowBkWvyT4idkuaBywHGoAlEbFK0kKgIyLagHmSzgZ2AW8Ds9PmZwALJe0G9gBXRcSWwXghZmZWmeptbKqlpSU6Ojp6r2hmQ67aPXvPgNl/kp6MiJZK2z1dgplZATjZm5kVgJO9mVkBONmbmRWAk72ZWQE42ZuZFYCTvZlZATjZm5kVgJO9mVkBONmbmRWAk72ZWQE42ZuZFYCTvZlZATjZm5kVgJO9mVkBONmbmRVAntsSmpntY6TdjrAIcvXsJU2XtEZSp6T5ZbZfJek5SU9L+o2kyZltN6bt1kg6r5rBm9nI0rqi9YOHVVevyV5SA7AIOB+YDFycTeap+yPiExFxCnA7cEfadjLJDcpPBKYD3033Z2ZmNZSnZz8V6IyIdRHxPrAUmJmtEBFbM6sHA903tp0JLI2InRHxMtCZ7s/MzGooz5j9eGBDZr0LOK20kqSrgeuAUcCZmbYrS9qO71ekZmbWb3l69ipTFvsVRCyKiOOAG4Cb+9JW0lxJHZI6Nm3alCMkMzPrizzJvguYkFlvBDb2UH8pcEFf2kbE4ohoiYiWcePG5QjJzMz6Ik+ybweaJTVJGkVywrUtW0FSc2b1L4G16XIbMEvSQZKagGbgtwMP28zM+qLXMfuI2C1pHrAcaACWRMQqSQuBjohoA+ZJOhvYBbwNzE7brpL0ILAa2A1cHRF7Bum1mJlZBbl+VBURy4BlJWULMsvX9tD2VuDW/gZoZmYD5+kSzMwKwMnezKwAnOzNzArAyd7MrACc7M3MCsDJ3sysAJzszcwKwMnezKwAfKcqM6tL2RuYtE5rrVjP8nHP3sysAJzszcwKwMnezKwAnOzNzArAyd7MrAB8NY6Z5ZK9OsaGH/fszcwKwMnezKwAciV7SdMlrZHUKWl+me3XSVot6VlJv5Z0TGbbHklPp4+20rZmZjb4eh2zl9QALALOAbqAdkltEbE6U+0poCUi3pP0FeB24Evptu0RcUqV4zYzsz7I07OfCnRGxLqIeB9YCszMVoiIxyLivXR1JdBY3TDNzGwg8iT78cCGzHpXWlbJHODhzPpoSR2SVkq6oFwDSXPTOh2bNm3KEZKZmfVFnksvVaYsylaULgNagM9miidGxEZJxwKPSnouIl7aZ2cRi4HFAC0tLWX3bWZm/ZenZ98FTMisNwIbSytJOhu4CZgRETu7yyNiY/rvOmAFMGUA8ZqZWT/kSfbtQLOkJkmjgFnAPlfVSJoC3E2S6N/KlB8h6aB0eSzwGSB7YtfMzGqg12GciNgtaR6wHGgAlkTEKkkLgY6IaAO+BRwC/FgSwKsRMQP4OHC3pL0kHyzfLLmKx8zMaiDXdAkRsQxYVlK2ILN8doV2/wZ8YiABmpnZwPkXtGZmBeBkb2ZWAE72ZmYF4GRvZlYATvZmZgXgZG9mVgBO9mZmBeBkb2Z1r3VFq2+LOEBO9mZmBeAbjptZRe5Njxzu2ZuZFYCTvZlZATjZm5kVgJO9mVkBONmbmRWAk72ZWQE42ZuZFUCuZC9puqQ1kjolzS+z/TpJqyU9K+nXko7JbJstaW36mF3N4M3MLJ9ek72kBmARcD4wGbhY0uSSak8BLRFxMvAQcHva9kjg68BpwFTg65KOqF74ZmaWR56e/VSgMyLWRcT7wFJgZrZCRDwWEe+lqyuBxnT5POCRiNgSEW8DjwDTqxO6mZnllWe6hPHAhsx6F0lPvZI5wMM9tB1f2kDSXGAuwMSJE3OEZNZHra19KzcbYfIke5Upi7IVpcuAFuCzfWkbEYuBxQAtLS1l922WmxO42X7yDON0ARMy643AxtJKks4GbgJmRMTOvrQ1M7PBladn3w40S2oCXgNmAZdkK0iaAtwNTI+ItzKblgN/nzkpey5w44CjNquWct8C/M3ARqBek31E7JY0jyRxNwBLImKVpIVAR0S0Ad8CDgF+LAng1YiYERFbJN1C8oEBsDAitgzKKzEzs4pyzWcfEcuAZSVlCzLLZ/fQdgmwpL8BmlXkHrhZbr55iZnto55vWJKNrXVaa8V6tj9Pl2BmVgDu2ZuV8klbG4HcszczKwAnezOzAvAwjg0PHkYxGxD37M3MCsDJ3sysAJzszcwKwGP2Znl4imQb5tyzNzMrACd7M7MCcLI3MysAJ3szswJwsjczKwBfjWP1x1e4mFWde/ZmZgWQq2cvaTrwbZLbEv5zRHyzZPsZwD8CJwOzIuKhzLY9wHPp6qsRMaMagZvVBU+HbMNEr8leUgOwCDgH6ALaJbVFxOpMtVeBK4Dry+xie0ScUoVYzcw+4LtW9U2env1UoDMi1gFIWgrMBD5I9hGxPt22dxBiNDOzAcozZj8e2JBZ70rL8hotqUPSSkkXlKsgaW5ap2PTpk192LWZmeWRJ9mrTFn04TkmRkQLcAnwj5KO229nEYsjoiUiWsaNG9eHXZuZWR55kn0XMCGz3ghszPsEEbEx/XcdsAKY0of4zMysCvKM2bcDzZKagNeAWSS99F5JOgJ4LyJ2ShoLfAa4vb/BmtngyJ7stJGp1559ROwG5gHLgReAByNilaSFkmYASDpVUhfwReBuSavS5h8HOiQ9AzwGfLPkKh4zM6uBXNfZR8QyYFlJ2YLMcjvJ8E5pu38DPjHAGG2k8vXoZjXjX9CamRWAk72ZWQF4IjSzavMtDK0OuWdvZlYATvZmZgXgZG9mVgBO9mZmBeBkb2ZWAE72ZmYF4GRvZlYATvZmZgXgH1VZbfgHRWZDysnerFZ8c3IbQk72ZgU2Uuax983He+cxezOzAnCyNzMrACd7M7MCyJXsJU2XtEZSp6T5ZbafIel3knZLuqhk22xJa9PH7GoFbmZm+fWa7CU1AIuA84HJwMWSJpdUexW4Ari/pO2RwNeB04CpwNfTm5CbmVkN5bkaZyrQGRHrACQtBWYCH9w4PCLWp9v2lrQ9D3gkIrak2x8BpgM/GnDkVr98OaFZ3ckzjDMe2JBZ70rL8sjVVtJcSR2SOjZt2pRz12ZmlleeZK8yZZFz/7naRsTiiGiJiJZx48bl3LWZmeWVJ9l3ARMy643Axpz7H0hbMzOrkjzJvh1oltQkaRQwC2jLuf/lwLmSjkhPzJ6blpmZWQ31eoI2InZLmkeSpBuAJRGxStJCoCMi2iSdCvwUOAL4gqRvRMSJEbFF0i0kHxgAC7tP1poZni/HaibX3DgRsQxYVlK2ILPcTjJEU67tEmDJAGI0M7MB8i9ozcwKwMnezEaU1hWtI2Y2z2ryFMdmBeNEWEzu2ZuZFYB79tZ/vmrEbNhwz97MrACc7M3MCsDJ3sysADxmb1ZvKp0L8TkSGwD37M3MCsDJ3sysAJzszcwKwMnezKwAfILW8vHJQRtmstNCtE5rrVivKJzszQrA8+GYk73ZcOEbndgAeMzezKwAciV7SdMlrZHUKWl+me0HSXog3f6EpElp+SRJ2yU9nT7uqm74ZmaWR6/DOJIagEXAOUAX0C6pLSJWZ6rNAd6OiOMlzQJuA76UbnspIk6pctxmZtYHeXr2U4HOiFgXEe8DS4GZJXVmAvemyw8BZ0lS9cI0M7OByHOCdjywIbPeBZxWqU5E7Jb0LjAm3dYk6SlgK3BzRPzv0ieQNBeYCzBx4sQ+vQCrMp/wMxuR8vTsy/XQI2ed14GJETEFuA64X9KH96sYsTgiWiKiZdy4cTlCMjOzvsjTs+8CJmTWG4GNFep0SToAOAzYEhEB7ASIiCclvQScAHQMNHAzs7z8A6t8Pft2oFlSk6RRwCygraROGzA7Xb4IeDQiQtK49AQvko4FmoF11QndzMzy6rVnn47BzwOWAw3AkohYJWkh0BERbcD3gR9I6gS2kHwgAJwBLJS0G9gDXBURWwbjhZgVUg8/tPKvZi0r1y9oI2IZsKykbEFmeQfwxTLtfgL8ZIAxmpnZAHm6hCLzlTdmheHpEszMCsDJ3sysADyMYzbSpMNz09av2Kd4xRXTah6K1Q/37M3MCsA9e7MRZkVJj972VdQfWDnZF4WvvDErNA/jmJkVgHv2ZgUx7Z4V+5X5pG1xONmbWWEVafzeyd5sBPBJWeuNk/1I4xOxZlaGk71ZgZUbxweP5Y9EvhrHzKwAlNxMqn60tLRER4dvZJWLh2wKrdbj9EXp7Q/XE7WSnoyIlkrb3bM3MysA9+yHC/fijfq76mak9/aHUy+/t559rhO0kqYD3ya5LeE/R8Q3S7YfBNwH/DmwGfhSRKxPt90IzCG5LeE1EbG8H6+jOJzUzerGSLoOv9eefXrD8P8LnAN0kdyA/OKIWJ2p81Xg5Ii4StIs4MKI+JKkycCPgKnAx4BfASdExJ5Kz1eonr0Tu/Wg3nrxAzWSvgXUY+KvRs9+KtAZEevSHS4FZgKrM3VmAq3p8kPAnZKUli+NiJ3Ay+kNyacCj/f1hdQlJ2urgpGW1CupdJlnOfX+wTAce/x5kv14YENmvQs4rVKdiNgt6V1gTFq+sqTt+NInkDQXmJuubpO0JkdcY4Hf56g3FOo5Nqjv+Bxb/9VzfH2L7d5/GbxIyuv3sfsG36hyKPvJG9sxPW3Mk+xVpqx07KdSnTxtiYjFwOIcsfz/J5Q6evrKMpTqOTao7/gcW//Vc3z1HBvUd3zVii3PpZddwITMeiOwsVIdSQcAhwFbcrY1M7NBlifZtwPNkpokjQJmAW0lddqA2enyRcCjkZz5bQNmSTpIUhPQDPy2OqGbmVlevQ7jpGPw84DlJJdeLomIVZIWAh0R0QZ8H/hBegJ2C8kHAmm9B0lO5u4Gru7pSpw+6tOwT43Vc2xQ3/E5tv6r5/jqOTao7/iqElvd/ajKzMyqz9MlmJkVgJO9mVkB1HWyl/RFSask7ZXUUrLtRkmdktZIOq9C+yZJT0haK+mB9ATzYMT5gKSn08d6SU9XqLde0nNpvZr9TFhSq6TXMjF+vkK96enx7JQ0v0axfUvSi5KelfRTSYdXqFezY9fbcUgvOHgg3f6EpEmDGU/Jc0+Q9JikF9L/G9eWqTNN0ruZv/eCGsbX499Jie+kx+5ZSZ+sUVx/mjkeT0vaKulvS+rU9LhJWiLpLUnPZ8qOlPRImrMekXREhbaz0zprJc0uV2c/EVG3D+DjwJ8CK4CWTPlk4BngIKAJeAloKNP+QWBWunwX8JUaxPxfgQUVtq0Hxg7BcWwFru+lTkN6HI8FRqXHd3INYjsXOCBdvg24bSiPXZ7jAHwVuCtdngU8UMO/5dHAJ9PlQ0mmMimNbxrw81q/z/L8nYDPAw+T/AbnU8ATQxBjA/AGcMxQHjfgDOCTwPOZstuB+eny/HL/H4AjgXXpv0eky0f09nx13bOPiBciotyvaT+YhiEiXga6p2H4QDpdw5kk0zcA3AtcMJjxps/5H0jmAxpuPpgWIyLeB7qnxRhUEfHLiNidrq4k+S3GUMpzHGaSvJ8geX+dlf7tB11EvB4Rv0uX/wC8QJlfpdexmcB9kVgJHC7p6BrHcBbwUkS8UuPn3UdE/CvJ1YtZ2fdWpZx1HvBIRGyJiLeBR4DpvT1fXSf7HpSbwqH0DT8GeCeTSMpO1VBlfwG8GRFrK2wP4JeSnkyniKileenX5iUVvhrmOaaD7UqSXl85tTp2eY7DPtODAN3Tg9RUOnw0BXiizOZPS3pG0sOSTqxhWL39nerhfTaLyh2yoTpu3Y6KiNch+WAHPlKmTr+O4ZDfg1bSr4CPltl0U0T8z0rNypTlncKhX3LGeTE99+o/ExEbJX0EeETSi+mn+4D1FB/wPeAWktd/C8lQ05WluyjTtirX5eY5dpJuIvktxg8r7GbQjl1puGXKBvW91R+SDgF+AvxtRGwt2fw7kiGKben5mZ+R/KCxFnr7Ow3psUvP280AbiyzeSiPW1/06xgOebKPiLP70SzPNAy/J/mKeEDa+xrQVA29xalkmoi/JpnTv9I+Nqb/viXppyRDBlVJWHmPo6R/An5eZtOgTW2R49jNBv4KOCvSQcky+xi0Y1eiL9ODdGnf6UFqQtKBJIn+hxHxP0q3Z5N/RCyT9F1JYyNi0CdJy/F3GuopVM4HfhcRb5ZuGMrjlvGmpKMj4vV0eOutMnW6SM4vdGskOa/Zo+E6jNPrNAxp0niMZPoGSKZzqPRNoRrOBl6MiK5yGyUdLOnQ7mWSE5PPl6tbbSVjohdWeN4802IMRmzTgRuAGRHxXoU6tTx2A5keZNCl5wa+D7wQEXdUqPPR7nMIkqaS/D/fXIPY8vyd2oDL06tyPgW82z1sUSMVv30P1XErkX1vVcpZy4FzJR2RDsmem5b1rFZnnvt5tvpCkk+xncCbwPLMtptIrppYA5yfKV8GfCxdPpbkQ6AT+DFw0CDGeg9wVUnZx4BlmVieSR+rSIYwanUcfwA8BzybvpmOLo0vXf88ydUdL9UqvvRvswF4On3cVRpbrY9dueMALCT5QAIYnb6fOtP317E1/Fv+e5Kv7M9mjtnngau633/AvPQ4PUNy0vv0GsVW9u9UEpuARemxfY7MVXY1iO9PSJL3YZmyITtuJB86rwO70jw3h+Tcz6+Btem/R6Z1W0juEtjd9sr0/dcJ/Mc8z+fpEszMCmC4DuOYmVkfONmbmRWAk72ZWQE42ZuZFYCTvZlZATjZm5kVgJO9mVkB/D/VaN1LAduNJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])\n",
    "\n",
    "def p_ygivenx(x, m1, m2, s1, s2):\n",
    "    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))\n",
    "\n",
    "def p_xgiveny(y, m1, m2, s1, s2):\n",
    "    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))\n",
    "\n",
    "N = 5000\n",
    "K = 20\n",
    "x_res = []\n",
    "y_res = []\n",
    "z_res = []\n",
    "m1 = 5\n",
    "m2 = -1\n",
    "s1 = 1\n",
    "s2 = 2\n",
    "\n",
    "rho = 0.5\n",
    "y = m2\n",
    "\n",
    "for i in range(N):\n",
    "    for j in range(K):\n",
    "        x = p_xgiveny(y, m1, m2, s1, s2)   #y给定得到x的采样\n",
    "        y = p_ygivenx(x, m1, m2, s1, s2)   #x给定得到y的采样\n",
    "        z = samplesource.pdf([x,y])\n",
    "        x_res.append(x)\n",
    "        y_res.append(y)\n",
    "        z_res.append(z)\n",
    "\n",
    "num_bins = 50\n",
    "plt.hist(x_res, num_bins,density=1, facecolor='green', alpha=0.5,label='x')\n",
    "plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5,label='y')\n",
    "plt.title('Histogram')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "OfldxGWecECT"
   },
   "source": [
    "----\n",
    "本章代码来源：https://github.com/hktxt/Learn-Statistical-Learning-Method\n",
    "\n",
    "本文代码更新地址：https://github.com/fengdu78/lihang-code\n",
    "\n",
    "中文注释制作：机器学习初学者公众号：ID:ai-start-com\n",
    "\n",
    "配置环境：python 3.5+\n",
    "\n",
    "代码全部测试通过。\n",
    "![gongzhong](../gongzhong.jpg)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "MCMC.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
